# Numerical simulation changes with different tspan

**URL:** <https://discourse.julialang.org/t/numerical-simulation-changes-with-different-tspan/102726>\
**Category:** Numerics\
**Tags:** differentialequation\
**Created:** [August 11, 2023, 5:43pm UTC](https://discourse.julialang.org/t/numerical-simulation-changes-with-different-tspan/102726 "2023-08-11T17:43:40Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Yoni\_Maltsman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yoni_maltsman/32/51557_2.png) [@Yoni\_Maltsman](https://discourse.julialang.org/u/Yoni_Maltsman)\
**Post date:** [August 11, 2023, 5:43pm UTC](https://discourse.julialang.org/t/numerical-simulation-changes-with-different-tspan/102726/1 "2023-08-11T17:43:40Z")

</div>

I am attempting to simulate a network of four neurons, where the `i`th neuron’s dynamics are described by the following differential equations:

\dot{V}\_m = I\_{syn} - (V\_m - \alpha\_f^-\tanh(V\_m - \delta\_f^-) + \alpha\_s^+\tanh(V\_s - \delta\_s^+) - \alpha\_s^-\tanh(V\_s - \delta\_s^-) + \alpha\_{us}^+\tanh(V\_{us} - \delta\_{us}^+))

\tau\_s \dot{V}\_s = V\_m - V\_s

\tau\_{us} \dot{V}\_{us} = V\_m - V\_{us}

where \tau\_{us} \>\> \tau\_{s} \>\> 1, and

I\_{syn} = I\_{app} + \sum\_{j \neq i} \frac{\alpha\_{ij}}{1+\exp(-\beta(V\_j^s-\delta\_{syn})}

The network is entirely homogeneous, where each neuron has the same parameters, and all of the weights between the neurons are the same.

I am using Lsoda and DifferentialEquations.JL to solve the system, with reltol = abstol = 1e-8. I find that the system either immediately converges to a fixed point or a limit cycle depending on the time span over which I integrate. Any suggestions as to why this is happening and how I might fix it would be appreciated.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [August 11, 2023, 5:48pm UTC](https://discourse.julialang.org/t/numerical-simulation-changes-with-different-tspan/102726/2 "2023-08-11T17:48:18Z")

</div>

Can you provide some Julia code?

---

<div class="post-metadata">

**Author:** ![Yoni\_Maltsman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yoni_maltsman/32/51557_2.png) [@Yoni\_Maltsman](https://discourse.julialang.org/u/Yoni_Maltsman)\
**Post date:** [August 11, 2023, 5:58pm UTC](https://discourse.julialang.org/t/numerical-simulation-changes-with-different-tspan/102726/3 "2023-08-11T17:58:17Z")

</div>

Sure. Here is the function for the ODE for the entire network:

```julia
function NNdynamics!(dV, V, p, t)
		@unpack neuronp, edgep, W = p
		@unpack β, δ = edgep
		for j in range(1, length = length(neuronp)) #loop through each neuron
			idx = (j-1)*3
			Vneuron = [V[idx+1], V[idx+2], V[idx+3]]
			dV[idx+1], dV[idx + 2], dV[idx + 3] = neuronDS(Vneuron, neuronp[j]) #internal dynamics
			for i in 1:length(W[j]) #again, loop through each neuron
				dV[idx+1] += W[i][j]*sigmoid(β*(V[((i-1)*3)+2] - δ)) #Isyn
			end
		end
	end

```

Here is the function for the internal dynamics of an individual neuron:

```julia
function neuronDS!(dV, V, p, t)
		@unpack αfm, αsp, αsm, δfm, δsp, δsm, αusp, δusp, Iₐ = p
		Vm, Vs, Vus = V
		dV[1] = Iₐ - ( Vm - αfm*tanh(Vm - δfm) + αsp*tanh(Vs - δsp) - αsm*tanh(Vs - δsm) + αusp*tanh(Vus - δusp))
		dV[2] = (1/τs)*(Vm - Vs)
		dV[3] = (1/τus)*(Vm - Vus)
		dV
	end

```

And here is my function for simulating the system. The thing that I am confused about is that changing T gives me qualitatively different results

```julia
function numsim(f, V0, par, T)
		prob = ODEProblem(f, V0, (0,T), par)
		sol = solve(prob, lsoda(), reltol = 1e-8, abstol = 1e-8)
	end

```

Let me know if there is anything else I should provide.

---

<div class="post-metadata">

**Author:** ![moble](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/moble/32/23535_2.png) [@moble](https://discourse.julialang.org/u/moble)\
**Post date:** [August 14, 2023, 3:28pm UTC](https://discourse.julialang.org/t/numerical-simulation-changes-with-different-tspan/102726/4 "2023-08-14T15:28:00Z")

</div>

I’m not an expert, but since nobody else is jumping in, I’ll make a few points.

The first thing that jumps out at me is that when you change `T`, the initial stepsize changes. It may be the case that when `T` is larger, you basically just step over an important dynamical feature, for example. I can’t find any documentation about how that initial stepsize is chosen automatically, and though you can pass the `dt` option to choose just that initial stepsize, it appears that `lsoda` doesn’t actually use that option.

That leads me to wonder why you’re using `lsoda`. In the first post, you pointed out that \tau\_{us} \gg \tau\_{s} \gg 1, which looks like you’re saying that the dynamics work on very different timescales, which suggests a stiff system. It is true that `lsoda` has long been a workhorse for that kind of problem, but it looks like there are typically better options nowadays — and certainly options that have better integration with the Julia interface. You could just use

```julia
sol = solve(prob, alg_hints = [:stiff], reltol = 1e-8, abstol = 1e-8)

```

to let the library choose an algorithm for you. Or if you want to get into the details, you should look at [this page](https://docs.sciml.ai/DiffEqDocs/stable/getting_started/#Choosing-a-Solver-Algorithm) (or find even more detail on [this page](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Recommended-Methods)). If I’m right in thinking that this is a stiff system, it looks like `RadauIIA5` might be the algorithm for you because you’re also asking for pretty high accuracy. Note that these pages seem to suggest `lsoda` unenthusiastically; its main use case seems to be when the number of variables is very large, but even then `QNDF`, `FBDF`, and `CVODE_BDF` all come earlier in the recommendations — but you say you only have 4 neurons in any case.

So I would investigate with different algorithms, the `dt` parameter, and the tolerances, to see if the fixed point or limit cycle behaviors are robust.

---

<div class="post-metadata">

**Author:** ![Yoni\_Maltsman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yoni_maltsman/32/51557_2.png) [@Yoni\_Maltsman](https://discourse.julialang.org/u/Yoni_Maltsman)\
**Post date:** [August 14, 2023, 5:54pm UTC](https://discourse.julialang.org/t/numerical-simulation-changes-with-different-tspan/102726/5 "2023-08-14T17:54:37Z")

</div>

Thank you very much! I was using lsoda in order to reproduce the results of a colleague who had run the simulations in Python with scipy’s odeint function, which uses the lsoda algorithm. I will try messing with the algorithm and parameters as you mentioned to see if I can still reproduce his results and get more robust behavior with the limit cycle and fixed points.