Summary: When I simulate the stochastic time development of a reaction network the simulation crashes due to concentrations becoming negative (works fine for deterministic solution).

I have this chemical reaction network:

```
2, ∅ --> X
3, X --> ∅
```

(original system more complicated, but I have simplified it while retaining the error)

I model it as:

```
function f(du,u,p,t)
du[1] = 2 - 3u[1]
end
function g(du,u,p,t)
du[1, 1] = sqrt(2)
du[1, 2] = -1 * sqrt(3u[1])
end
```

when I simulate the system using the deterministic approach everything is fine:

```
u0 = [50.0]
tspan = (0.0,4.0)
prob_det = ODEProblem(f,u0,tspan)
sol_det = solve(prob_det)
```

however when I try to solve the stochastic version it fails:

```
u0 = [50.0]
tspan = (0.0,4.0)
dt = 1//2^(12)
prob_stoch = SDEProblem(f,g,u0,tspan,noise_rate_prototype=zeros(1,2))
sol_stoch = solve(prob_stoch,EM(),dt=dt)
```

the last line yields a `DomainError()`

which can be traced to `du[1, 2] = -1 * sqrt(3u[1])`

(where the value of `u[1]`

is negative).

I’ve spent some time trying to figure out how this happens, and I have no idea. I’ve modified the functions to print values:

```
function f(du,u,p,t)
println(u[1], " ",(2 - 3u[1]))
...
end
function g(du,u,p,t)
println(sqrt(2), " ",-1 * sqrt(3u[1]),"\n")
...
end
```

this yields:

```
...
0.16617661194477099 1.5014701641656871
1.4142135623730951 -0.7060664528458444
0.18415292624496235 1.447541221265113
1.4142135623730951 -0.7432757084251355
0.18156840309376518 1.4552947907187046
1.4142135623730951 -0.7380414685377072
0.2004141078247585 1.3987576765257246
1.4142135623730951 -0.7753981709252836
0.164105752037157 1.507682743888529
1.4142135623730951 -0.7016532306712989
0.1469872071015845 1.5590383786952464
1.4142135623730951 -0.6640494117945994
0.15114573738914047 1.5465627878325785
1.4142135623730951 -0.6733774663347604
0.1506813504673137 1.5479559485980587
1.4142135623730951 -0.672342213015025
0.11613339743580103 1.6515998076925968
1.4142135623730951 -0.5902543454371201
0.05865690190122183 1.8240292942963345
1.4142135623730951 -0.4194886240455937
0.061455047206756666 1.81563485837973
1.4142135623730951 -0.42937762123831047
0.03440580432279003 1.89678258703163
1.4142135623730951 -0.3212746690425034
0.01917417408706192 1.9424774777388143
1.4142135623730951 -0.23983853372881048
0.028964799664317847 1.9131056010070464
1.4142135623730951 -0.2947785592490633
0.026156451214798624 1.921530646355604
1.4142135623730951 -0.28012381841677775
-0.005493955043100477 2.0164818651293013
```

Here both the deterministic and stochastic part seems overwhelmingly positive (which makes sense as the degradation of X should be low as the concentration of X approaches 0). Yet the value of X is still being steadily reduced.

(I have been running this simulation a lot, and it always yields the same error after a while)

There is surely something here which I do not understand, I would be happy if someone could help me point out what that is.