DifferentialEquations: Gaussian error term in one of the equations makes it extremely slow (and error prone)


For my thesis, I’m trying to implement this model in Julia, using the DifferentialEquations.jl package. I’ve made a rule based on a set of coupled differential equations to simulate some data. I have no previous knowledge of the subject, so bear in mind that I might have to awkwardly stumble over words to explain my problem:

function threePlusOneDimensions!(du, u, p, t)

    gauss = 0.2
    Smax, τx, τy, τz, τf, α, β, P, S, Rₛ, L, λᵦ, λₛ = p

    du[1] = (Smax /(1 + exp((Rₛ - u[2]) / λₛ)) - u[1]) / τx
    du[2] = (P / (1 + exp((Rₛ - u[2]) / λₛ)) + u[4] * L - u[1] * u[2] - u[3]) / τy
    du[3] = (S * (α * u[1] + β * u[2]) - u[3])  * gauss / τz
    du[4] = (u[2] - λᵦ * u[4])/ τf
    return nothing

This is the rule I’m specifying. The gauss variable should “integrate the equation as external noise ζ(t), set between –1 and 1 with Gaussian distribution.” Without the drawn value, it is very fast. As soon as I add the external noise term, it slows down tremendously. How do I implement this in an idiomatic manner? I copied the full text from the paper explaining the third equation below.

Equation 3: Modeling of perceived environment

The third equation refers to the environment (or external world) perceived by a patient, modeled by the variable z (Equation 3), which evolves with a time constant τz. It depends on the overal sensitivity level S, and the joint effects of symptoms x and the potentiation y respectively pondered by factor α and β. The factors α and β may be positive or negative depending on the type of psychiatric disease considered. The perceived environment integrates the equation as external noise ζ(t), set between –1 and 1 with Gaussian distribution. The release occurs with an exponential decay (−z).

This part of the paper is also relevant:

Thirdly, the rate of the noise ζ(t) is chosen at 0.01, meaning that the perceived environment variable z changes every 0.01 days (noise will be generated every 14.4 min). It is a compromise between the duration of variability of the symptoms of psychiatric disorders and their environment (i.e., considering a psychological state change every 14.4 min). In other words, the model provides a smoothness of 14.4 min, i.e., informs about potential changes in its variables approximately every quarter of an hour.

How do I do this idiomatically, so that the package properly recognizes that there is a random term? I’ve tried doing it with rand() or through drawing from a distribution, but I cannot do so according to the specifications in the paper without it slowing down tremendously. I wanted to solve it myself, and also tried many other things, asked my supervisors for help as well, but I haven’t been able to proceed for some time now.

You should implement this as a discrete callback which changes a parameter according this this schedule. I recommend using a PeriodicCallback.


Thank you! This is exactly the kind of clear hint I was looking for!