ModelingToolkit and Noise Processes for uniform noise definition

I tried another approach with symbolically registering rand() function and “artificially” adding noise equations to the ODESystem to create SDESystem, and seems this is a solution for my problem:

using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations
using Plots
using Interpolations

# Sigmoid function
function sigmoid(V, e0, r, θ)
    return 2 * e0 / (1 + exp(r * (θ - V)))
end

# Define a time-varying random function
function random_input(t)
    return 320.0 * rand()
end
@register_symbolic (random_input)(t)

params = @parameters A0=3.25 a0=100.0 A1=3.25 a1=100.0 A2=22.0 a2=50.0 C1=135.0 C2=0.8 * C1 C3=0.25 * C1 C4=0.25 * C1 e0=5.0 r=0.56 v0=6.0
vars = @variables (u(t))[1:6]

eqs = [D(u[1]) ~ u[4],
       D(u[2]) ~ u[5],
       D(u[3]) ~ u[6],
       D(u[4]) ~ A0 * a0 * sigmoid(u[2] - u[3], e0, r, v0) - 2 * a0 * u[4] - a0^2 * u[1],
       D(u[5]) ~ A1 * a1 * (random_input(t) + C2 * sigmoid(C1 * u[1], e0, r, v0)) - 2 * a1 * u[5] - a1^2 * u[2],
       D(u[6]) ~ A2 * a2 * (C4 * sigmoid(C3 * u[1], e0, r, v0)) - 2 * a2 * u[6] - a2^2 * u[3]
       ]

noiseeqs = [0., 0., 0., 0., 1., 0.]

# Parameters
params_map = Dict(A0 => 3.25, a0 => 100.0, A1 => 3.25, a1 => 100.0, 
                  A2 => 22.0, a2 => 50.0, C1 => 135.0, C2 => 0.8 * 135.0, 
                  C3 => 0.25 * 135.0, C4 => 0.25 * 135.0, e0 => 5.0, r => 0.56, v0 => 6.0)

inits = [0.39057, -204.28131, -3.88133, -0.11387, -7.68439, -253.32597]
tspan = (0.0, 10.0)
dts = 1e-2
dt = 1e-4

sys = structural_simplify(SDESystem(eqs, noiseeqs, t, vars, params; name=Symbol("JR")))
prob = SDEProblem(sys, inits, tspan, params_map) 
sol = solve(prob, EM(), dt=dt, saveat=dts)
plot(sol[1, :])