I am relatively new to Mathematica and I need some help regarding how to create a random number generator to be used as an independent variable for a model my group has been working on. Following is our current base model:

 sol[fmaxnW_?NumericQ, fminnW_?NumericQ, fmaxpW_?NumericQ, 
  fmaxnR_?NumericQ, fminnR_?NumericQ, micW_?NumericQ, micR_?NumericQ, 
  kW_?NumericQ, kR_?NumericQ , K_?NumericQ, a_?NumericQ, b_?NumericQ, 
  mu_?NumericQ, Amax_?NumericQ, w_?NumericQ, pinit_?NumericQ, 
  tinit_?NumericQ, propR_?NumericQ] :=
  s = NDSolve[
    {A'[t] == -w*A[t],

 nW'[t] == (fmaxnW - ((fmaxnW - fminnW)*(A[t]/micW)^
          kW)/((A[t]/micW)^kW - fminnW/fmaxnW))*
    Which[A[t] < micW, (1 - 1/K (nW[t] + pW[t] + nR[t])), True, 
     1]*nW[t]  - a*1/K (nW[t] + pW[t] + nR[t])*nW[t]  +  b*pW[t] -
    w*nW[t] - 
   mu*(fmaxnW - ((fmaxnW - fminnW)*(A[t]/micW)^kW)/((A[t]/micW)^
          kW - fminnW/fmaxnW))*nW[t],

 pW'[t] == 
  fmaxpW*pW[t]  +  a*1/K (nW[t] + pW[t] + nR[t])*nW[t]  -  
   b*pW[t] - w*pW[t] - mu*fmaxpW*pW[t],

 nR'[t] == (fmaxnR - ((fmaxnR - fminnR)*(A[t]/micR)^
          kR)/((A[t]/micR)^kR - fminnR/fmaxnR))*
    Which[A[t] < micR, (1 - 1/K (nW[t] + pW[t] + nR[t])), True, 
     1]*nR[t] - w*nR[t] + 
   mu*(fmaxnW - ((fmaxnW - fminnW)*(A[t]/micW)^kW)/((A[t]/micW)^
          kW - fminnW/fmaxnW))*nW[t] + mu*fmaxpW*pW[t],

 (* initial conditions *)
 A[0] == 0, WhenEvent[Mod[t, 24] == 0, A[t] -> Amax],
 pW[0] == pinit*(1 - propR), 
 nW[0] == (tinit - pinit)*(1 - propR),  nR[0] == tinit*propR
    {A, nW, pW, nR}, {t, 0, 240}];

sol[1, -6, 0.00001, 0.9, -6, 0.01, 0.05, 1, 1, 
 5*10^9, 0.0001, 0.14, mu, 0.02, 0.231, 1, 2*10^5, 0]

Where nW represents normal cells, pW represents persister cells, and nR represents resistant cells. The model is currently a deterministic approach to switching rate and mutation of resistance in bacteria populations and we wish to make mutation rate ('mu' in this model) a stocastic variable. Specifically, we want 'mu' to generate a new random number every time step (t) in order to represent random probability in mutation rate (in the range of: 6*10^(-8), 10*10^(-8)). How can we make mu as a random number generator and independent variable for this model? Thanks in advance!

  • $\begingroup$ What do you mean, "every time step"? NDSove is unlikely to use a fixed time step. Do you mean that you want mu to randomly change every once in a while (like, say, every $\Delta t = 0.1$ or something)? $\endgroup$
    – march
    Jan 4, 2016 at 2:34
  • $\begingroup$ @march Yes, that's what I was meant to say. $\endgroup$
    – WDemoen
    Jan 4, 2016 at 3:31
  • 2
    $\begingroup$ I don't really understand why you would want to change a parameter randomly in the middle of an integration. If you are looking for the sensitivity of your calculation to the probability, or what the spread is over a range of plausible values for the probability, then you should use a Monte Carlo approach. You would run the calculation many times, say a few thousand, where each time you use one fixed value of the probability that you selected randomly from the plausible distribution. Then you can look at the spread of results from the calculation $\endgroup$
    – Mark Adler
    Jan 4, 2016 at 6:09
  • $\begingroup$ Are you looking to do Ito/Stochastic Calculus? Are you looking to do simulation with something like the Euler-Maruyama method? en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method $\endgroup$
    – Searke
    Jan 4, 2016 at 15:14

1 Answer 1


Your code is rather long, so here is a simple example of how to implement a parameter that randomly changes every $\Delta t$.

dt = 0.2;
sol = NDSolve[{
     x'[t] == mu[t] x[t]
     , x[0] == 1
     , mu[0] == -1
     , WhenEvent[Mod[t, dt] == 0., mu[t] -> RandomReal[{-1, 1}]]
    , {x, mu}
    , {t, 0, 10}
    , DiscreteVariables -> {mu}

The trick is to promote the quantity mu to a DiscreteVariable and update its value every dt using WhenEvent inside NDSolve. You would change the quantities inside RandomReal[{-1, 1}] to whatever range you need for the random values of mu. Running the above code once and plotting yields

Plot[x[t] /. First@sol, {t, 0, 10}]

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.