improve SDE simulation desc

Author: TorkelE

docs/src/model_simulation/simulation_introduction.md

@@ -187,23 +187,34 @@ sol = solve(sprob, STrapezoid())
 sol = solve(sprob, STrapezoid(); seed = 123) # hide
 plot(sol)
 ```
-we can see that while this simulation (unlike the ODE ones) exhibits some fluctuations. Next, we reduce species amounts (using [`remake`](@ref ref)), thereby also increasing the relative amount of noise, we encounter a problem when the model is simulated:
+we can see that while this simulation (unlike the ODE ones) exhibits some fluctuations. 
+
+!!! note
+    Unlike for ODE and jump simulations, there are no good heuristics for automatically selecting suitable SDE solvers. Hence, for SDE simulations a solver must be provided. `STrapezoid` will work for a large number of cases. When this is not the case, however, please check the list of [available SDE solvers](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/) for a suitable alternative (making sure to select one compatible with non-diagonal noise and the [Ito interpretation]https://en.wikipedia.org/wiki/It%C3%B4_calculus).
+
+### [Common SDE simulation pitfalls](@id simulation_intro_SDEs_pitfalls)
+Next, let us reduce species amounts (using [`remake`](@ref ref)), thereby also increasing the relative amount of noise, we encounter a problem when the model is simulated:
 ```@example simulation_intro_sde
 sprob = remake(sprob; u0 = [:X1 => 0.33, :X2 => 0.66])
 sol = solve(sprob, STrapezoid())
 sol = solve(sprob, STrapezoid(); seed = 1234567) # hide
 plot(sol)
 ```
-Here, we receive a warning that the simulation was terminated before the designated final time point, also suggesting that an instability might have been encountered. What has happened here is that the stochastic fluctuations have grown too large for the adaptive time-stepper to handle. This can be remedied by disabling adaptive time-stepping and setting a low `dt`:
+Here, we receive a warning that the simulation was aborted. In the plot, we also see that it is incomplete. In this case we also note that species concentrations are very low (and sometimes, due to the relatively high amount of noise, even negative). This, combined with the early termination, suggests that we are simulating our model for too low species concentration for the assumptions of the CLE to hold. Instead, [jump simulations](@ref simulation_intro_jumps) should be used.
+
+Next, let us consider a simulation for another parameter set:
 ```@example simulation_intro_sde
-sol = solve(sprob, STrapezoid(); adaptive = false, dt = 0.001)
-sol = solve(sprob, STrapezoid(); seed = 1234567, adaptive = false, dt = 0.001) # hide
+sprob = remake(sprob; u0 = [:X1 => 100.0, :X2 => 200.0], p = [:k1 => 200.0, :k2 => 500.0])
+sol = solve(sprob, STrapezoid())
+sol = solve(sprob, STrapezoid(); seed = 12345) # hide
+plot(sol)
+```
+Again, the simulation is aborted. This time, however, species concentrations are relatively large, so the CLE might still hold. What has happened this time is that the accuracy of the simulations has not reached its desired threshold. This can be deal with [by reducing simulation tolerances](@ref ref):
+```@example simulation_intro_sde
+sol = solve(sprob, STrapezoid(), abstol = 1e-1, reltol = 1e-1)
+sol = solve(sprob, STrapezoid(); seed = 12345, abstol = 1e-1, reltol = 1e-1) # hide
 plot(sol)
 ```
-However, SDE solver instability (combined with the fact that the fluctuations push the solution to negative values!) suggests that species concentrations (for this combination of model, initial conditions, and parameter values) are too low for the CLE's approximating assumptions to hold. Instead, [jump simulations](@ref simulation_intro_jumps) should be used.
-
-!!! note
-    Unlike for ODE and jump simulations, there are no good heuristics for automatically selecting suitable SDE solvers. Hence, for SDE simulations a solver must be provided. `STrapezoid` will work for a large number of cases. When this is not the case, however, please check the list of [available SDE solvers](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/) for a suitable alternative (making sure to select one compatible with non-diagonal noise and the [Ito interpretation]https://en.wikipedia.org/wiki/It%C3%B4_calculus).
 
 ### [Scaling the noise in the chemical Langevin equation](@id simulation_intro_SDEs_noise_saling)
 When using the CLE to generate SDEs from a CRN, it can sometimes be desirable to scale the magnitude of the noise. This can be done by introducing a *noise scaling term*, with each noise term generated by the CLE being multiplied with this term. A noise scaling term can be set using the `@default_noise_scaling` option:
@@ -213,11 +224,11 @@ two_state_model = @reaction_network begin
   (k1,k2), X1 <--> X2
 end
 ```
-Here, we set the noise scaling term to `0.1`, reducing the noise with a factor $10$ (noise scaling terms $>1.0$ increase the noise, while terms $<1.0$ reduce the noise). If we simulate the modified model using the low-concentration inputs used previously, we see that the noise has been reduced (in fact by so much that the model can now be simulated without issues):
+Here, we set the noise scaling term to `0.1`, reducing the noise with a factor $10$ (noise scaling terms $>1.0$ increase the noise, while terms $<1.0$ reduce the noise). If we re-simulate the model using the low-concentration settings used previously, we see that the noise has been reduced (in fact by so much that the model can now be simulated without issues):
 ```@example simulation_intro_sde
-u0 = [:X1 => 0.33, :X2 => 0.66]
+u0 = [:X1 => 100.0, :X2 => 200.0]
 tspan = (0.0, 1.0)
-ps = [:k1 => 2.0, :k2 => 5.0]
+ps = [:k1 => 200.0, :k2 => 500.0]
 sprob = SDEProblem(two_state_model, u0, tspan, ps)
 sol = solve(sprob, STrapezoid())
 sol = solve(sprob, STrapezoid(); seed = 123) # hide