Simulation of the stochastic system

Question

I have the stochastic system which consists of 4 nonlinear equations. White Gaussian noise is used in the third equation only. Nevertheless, the whole system is stochastic.

Some problems arise when I try to compute random trajectories:

General::nomem: The current computation was aborted because there was insufficient memory available to complete the computation.

Throw::sysexc: Uncaught SystemException returned to top level. Can be caught with Catch[…, _SystemException].

There is plenty of free memory, but for some reasons MATHEMATICA considers my system as an uncomputable one. Is there a way to overcome it? here is my code:

(*parameters*)
ρ = 0.2; νmax = 100; Vth = 25; kV = 20; gKLeakOvergl = 0.5; \
GsynOvergl = 5; Kbath = 8.5; τK = 100;  γ = 10; δK \
= 0.02; τNa = 20;  δNa = 0.03; τm = 0.01; τD = \
2; σOvergl = 25; δxD = 0.01;
  
(*additional functions*)
    Ipump[K_, Na_] = ρ/((1 + Exp[3.5 - K]) (1 + Exp[(25 - Na)/3]));
    f[x_?NumericQ] = 
      Dot @@ MapAt[Boole, 
        Internal`FromPiecewise[Piecewise[{{0, x < 0}, {x, x >= 0}}]], 1];
    ν[V_] := νmax f[2/(1 + Exp[-2 (V - Vth)/kV]) - 1];
    Vk[K_] = 26.6 Log[K/130];
    uOvergL[K_, V_, xD_] = 
      gKLeakOvergl (Vk[K] - Vk[K0]) + GsynOvergl ν[V] (xD - 0.5); 

(*initial conditions*)
K0 = 3;  Na0 = 10; V0 = Vk[K0]; xD0 = 1;

(*differential equations*)
rhs1[K_, Na_, V_, xD_] := (Kbath - K)/τK - 
  2 γ Ipump[K, Na] + δK ν[V]; 
rhs2[K_, Na_, V_, xD_] := (Na0 - Na)/τNa - 
  3 Ipump[K, Na] + δNa ν[V];
rhs3[K_, Na_, V_, xD_] := (-V + uOvergL[K, V, xD])/τm;
rhs4[K_, Na_, V_, xD_] := (1 - xD)/τD - δxD xD ν[V];

proc = ItoProcess[{\[DifferentialD]K[t] == 
    rhs1[K[t], Na[t], V[t], 
      xD[t]] \[DifferentialD]t, \[DifferentialD]Na[t] == 
    rhs2[K[t], Na[t], V[t], 
      xD[t]] \[DifferentialD]t, \[DifferentialD]V[t] == 
    rhs3[K[t], Na[t], V[t], 
       xD[t]] \[DifferentialD]t + σOvergl/τm \
\[DifferentialD]w[t], \[DifferentialD]xD[t] == 
    rhs4[K[t], Na[t], V[t], xD[t]] \[DifferentialD]t}, {K[t], Na[t], 
   V[t], xD[t]}, {{K, Na, V, xD}, {K0, Na0, V0, xD0}}, 
  t, {w \[Distributed] WienerProcess[]}]

solproc = 
  RandomFunction[proc, {0., 100., 0.1}, Method -> "EulerMaruyama"];

Answer 1

First, I changed K to k, since K is a built-in symbol. Then I ran the process for a shorter time and found that V[t] got huge:

solproc = RandomFunction[proc, {0., 1.2, 0.1}, Method -> "EulerMaruyama"];
ListLinePlot[solproc, PlotRange -> All, PlotLegends -> {k, Na, V, xD}]

I then tried solving the deterministic part with NDSolve, which worked fine:

sol = NDSolve[{
  k'[t] == rhs1[k[t], Na[t], V[t], xD[t]], 
  Na'[t] == rhs2[k[t], Na[t], V[t], xD[t]], 
  V'[t] == rhs3[k[t], Na[t], V[t], xD[t]], 
  xD'[t] == rhs4[k[t], Na[t], V[t], xD[t]], k[0] == K0, 
  Na[0] == Na0, V[0] == V0, xD[0] == xD0}, {k, Na, V, xD}, {t, 0, 10}];

Plot[Evaluate[{k[t], Na[t], V[t], dX[t]} /. sol], {t, 0, 10}]

Going back to the stochastic model, decreasing the step size seems to fix the problem.

solproc = RandomFunction[proc, {0., 100, 0.01}, Method -> "EulerMaruyama"];
ListLinePlot[solproc, PlotRange -> All, PlotLegends -> {k, Na, V, xD}]

No apparent improvements with a smaller step size:

solproc = RandomFunction[proc, {0., 100, 0.001}, Method -> "EulerMaruyama"];
ListLinePlot[solproc, PlotRange -> All, PlotLegends -> {k, Na, V, xD}]