Differential Equations in Julia

  1. Lotka-Volterra system
  2. Izhikevich Neuron
    1. Simulating spikes
    2. Changing parameter values and initial conditions
    3. Adding external currents

Differential Equations in Julia

Lotka-Volterra system

# Import the packages we need
using ModelingToolkit ## for symbolic system definition
using ModelingToolkit: t_nounits as t, D_nounits as D ## generic symbolic time variable and derivative operator
using OrdinaryDiffEq ## for ODE solvers
using CairoMakie ## for plotting

# Time-dependent variables.
# The assigned values are default values that will be used as initial conditions if initial conditions are not explicitly provided.
@variables x(t)=1 y(t)=1
# Parameters with assigned values
@parameters a=1.5 b=1.0 c=3.0 d=1.0

# Model equations
eqs = [D(x) ~ a * x - b * x * y,
       D(y) ~ -c * y + d * x * y]

# Model system, defined by the equations, the independent time variable, the time-dependent variables and the parameters
@named sys = ODESystem(eqs, t, [x, y], [a, b, c, d])
# Simplified system equations. Necessary step before solving!
simpsys = structural_simplify(sys)

# Simulation timespan
tspan = (0.0, 10.0)

# Initial conditions. A vector of pairs of the form `[state => initial_condition, ...]`
u0 = [x => 5, y => 2]

# Problem to be solved
prob = ODEProblem(simpsys, u0, tspan)
# Solution of the problem using the Tsit5 solver
sol = solve(prob, Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 31-element Vector{Float64}:
  0.0
  0.08117449098753833
  0.2033472146318696
  0.34481923758030086
  0.5133453442586928
  0.6997205942636654
  0.9518449632884769
  1.2078720930549587
  1.5117501225731436
  1.8960367593015124
  2.3419689205858836
  2.8001611761311196
  3.154617596349518
  3.446603431988819
  3.816591344004686
  4.104070083015805
  4.50104637647137
  4.897856816920354
  5.496187010603642
  6.017436474469348
  6.362775057001467
  6.659006719074707
  6.972834263097586
  7.347100606008145
  7.7074811126632445
  8.176455290373381
  8.672007583256836
  9.165339638990314
  9.559050370550041
  9.893470554340231
 10.0
u: 31-element Vector{Vector{Float64}}:
 [5.0, 2.0]
 [4.737486951537538, 2.329442325406833]
 [4.159128176496565, 2.7852722096289324]
 [3.378232286801216, 3.105606216306078]
 [2.5672395472069955, 3.0811885263548477]
 [1.9773834953757778, 2.676854050573453]
 [1.609406579313759, 1.9582843875373555]
 [1.552749513812718, 1.3546873132985058]
 [1.7519768926099824, 0.8937306633628304]
 [2.3546312250537205, 0.61337657523862]
 [3.5563657586661104, 0.5897128591272462]
 [4.981843895194956, 1.0731410601672249]
 [4.818950256230912, 2.2392637119053838]
 [3.343661157987135, 3.1118008426116455]
 [1.9313561988828034, 2.6139960476246533]
 [1.5774359332463672, 1.8066110816035383]
 [1.6566423633822205, 1.02648877617623]
 [2.174120423530587, 0.6583968490231039]
 [3.750824823666615, 0.6146437130025499]
 [5.148338613885276, 1.3873773408422705]
 [4.315048391711306, 2.6760882747903865]
 [2.757183955074088, 3.1223343036598226]
 [1.814291633172061, 2.437547442159294]
 [1.5487493568378434, 1.4588265428851046]
 [1.761324481429152, 0.8887489752549641]
 [2.557619101510861, 0.5854753413658367]
 [4.022955232810912, 0.6610024111093685]
 [5.1545522666136785, 1.5571104452908116]
 [3.7679793462689086, 2.9724967230117327]
 [2.227954470813059, 2.900435675873252]
 [1.946115942732129, 2.6293333027601915]

The solution object contains the values of every variable of the system (x and y) for every simulated timestep. One can easily access the values of a specific variable using its symbolic name

sol[x] ## or similarly sol[y]
31-element Vector{Float64}:
 5.0
 4.737486951537538
 4.159128176496565
 3.378232286801216
 2.5672395472069955
 1.9773834953757778
 1.609406579313759
 1.552749513812718
 1.7519768926099824
 2.3546312250537205
 3.5563657586661104
 4.981843895194956
 4.818950256230912
 3.343661157987135
 1.9313561988828034
 1.5774359332463672
 1.6566423633822205
 2.174120423530587
 3.750824823666615
 5.148338613885276
 4.315048391711306
 2.757183955074088
 1.814291633172061
 1.5487493568378434
 1.761324481429152
 2.557619101510861
 4.022955232810912
 5.1545522666136785
 3.7679793462689086
 2.227954470813059
 1.946115942732129

Finally we can plot the timeseries of both variables of the solution in a line plot using

fig = plot(sol)
# display the figure
fig

We will shortly see more plot types and options.

Izhikevich Neuron

Simulating spikes

Let's now move on to a two-dimensional model of a neuron. This is a popular model, known as the Izhikevich neuron [1]. Even though the model appears simple, it can exhibit many different spiking patterns (e.g. regular and fast spiking, bursting, chattering) by changing its parameters. The Izhikevich neuron is a similar system of ODEs like the Lotka-Volterra example above. One notable difference however is the spiking mechanism. Spiking in the Izhikevich neuron needs to be implemented "manually". That is we need to detect when the voltage variable crosses a spiking threshold and every time this event happens we need to reset the neuron's voltage to a more polarized value and potentially alter other variables too.

@variables v(t)=-65 u(t)=-13
# The following parameter values correspond to chattering spiking.
@parameters a=0.02 b=0.2 c=-50 d=2 I=10

eqs = [D(v) ~ 0.04 * v ^ 2 + 5 * v + 140 - u + I,
        D(u) ~ a * (b * v - u)]
2-element Vector{Symbolics.Equation}:
 Differential(t)(v(t)) ~ 140 + I - u(t) + 5v(t) + 0.04(v(t)^2)
 Differential(t)(u(t)) ~ a*(-u(t) + b*v(t))

Now we need to define the callback to simulate spikes. Spiking callbacks belong to the family of discrete callbacks and are defined by two parts :

The spiking event for this model looks like this

event = (v > 30.0) => [u ~ u + d, v ~ c]
v(t) > 30.0 => Symbolics.Equation[u(t) ~ d + u(t), v(t) ~ c]

Now we can move on to define the system including the spiking event, create a problem and solve it.

# Model system and simplification.
@named izh_system = ODESystem(eqs, t, [v, u], [a, b, c, d, I]; discrete_events = event)
izh_simple = structural_simplify(izh_system)

# Initial conditions.
u0 = [v => -65.0, u => -13.0]

# Simulation timespan
tspan = (0.0, 400.0)

izh_prob = ODEProblem(izh_simple, u0, tspan)

izh_sol = solve(izh_prob)

fig = plot(izh_sol)
fig

or if we want to plot just the voltage timeseries with its spiking pattern

fig = plot(izh_sol; idxs=[v])
fig

Changing parameter values and initial conditions

After defining and simulating a system we might want to run another simulation by changing either or both of the parameter values and the initial conditions. The Izhikevich neuron that we simulated above has multiple spiking regimes that correspond to different parameter value combinations. If we have our ODEProblem already defined as the izh_prob variable above, we can remake it by keeping everything the same except for the fields that we choose to alter. Here we will change only the model parameters to change from chattering to fast spiking.

izh_prob = remake(izh_prob; p = [a => 0.1, b => 0.2, c => -65, d => 2])
# or we can change the initial conditions along with the parameters as
izh_prob = remake(izh_prob; p = [a => 0.1, b => 0.2, c => -65, d => 2], u0 = [v => -70, u => -10])

izh_sol = solve(izh_prob)

fig = plot(izh_sol; idxs=[v])
fig

Notice how the spiking pattern has changed compared to the previous simulation.

Adding external currents

Until now we have been simulating the Izhikevich neuron by injecting it with a constant external DC current I=10. We'll now see how we can expand the I input to a dynamic current, which is more realistic (currents do not remain constant in the brain for too long). Since I will change dynamically in time, it will be a time-dependent variable of the system and not a constant parameter.

@variables v(t)=-65 u(t)=-13 I(t)
# The following parameter values correspond to regular spiking.
@parameters a=0.02 b=0.2 c=-65 d=8

eqs = [D(v) ~ 0.04 * v ^ 2 + 5 * v + 140 - u + I,
        D(u) ~ a * (b * v - u),
        I ~ 10*sin(0.5*t)]

event = (v > 30.0) => [u ~ u + d, v ~ c]

# Notice how `I` was moved from the parameter list to the variable list in the following call.
@named izh_system = ODESystem(eqs, t, [v, u, I], [a, b, c, d]; discrete_events = event)
izh_simple = structural_simplify(izh_system)
Model izh_system:
Equations (2):
  2 standard: see equations(izh_system)
Unknowns (2): see unknowns(izh_system)
  v(t) [defaults to -65]
  u(t) [defaults to -13]
Parameters (4): see parameters(izh_system)
  a [defaults to 0.02]
  d [defaults to 8]
  b [defaults to 0.2]
  c [defaults to -65]
Observed (1): see observed(izh_system)

Let's display the equations of the original and the simplified system to see the effect of structural_simplify.

equations(izh_system)
3-element Vector{Symbolics.Equation}:
 Differential(t)(v(t)) ~ 140 - u(t) + 5v(t) + I(t) + 0.04(v(t)^2)
 Differential(t)(u(t)) ~ a*(-u(t) + b*v(t))
 I(t) ~ 10sin(0.5t)

The original system contains the algebraic equation for the external current I.

equations(izh_simple)
2-element Vector{Symbolics.Equation}:
 Differential(t)(v(t)) ~ 140.0 - u(t) + 5.0v(t) + I(t) + 0.04(v(t)^2)
 Differential(t)(u(t)) ~ a*(-u(t) + b*v(t))

However the I equation has been simplified away in izh_simple, since it was substituted into the differential equation D(v). This is an important functionality of structural_simplify as it has turned a system of Differential Algebraic Equations (DAE) to one of Ordinary Differential Equations (ODE). ODEs can be integrated (solved) much more efficiently compared to DAEs.

tspan = (0.0, 400.0)

# We do not provide initial conditions, so the default values for `v` and `u` from above will be used.
izh_prob = ODEProblem(izh_simple, [], tspan)

izh_sol = solve(izh_prob)

fig = plot(izh_sol; idxs=[v, I])
fig

Notice how the external current is slowly being accumulated in the neuron's potential v until the eventual spike and reset.

References

CC BY-SA 4.0 Neuroblox Inc. Last modified: December 09, 2024. Website built with Franklin.jl and the Julia programming language.