Implicit Nonlinear Dynamics : Autoregulation

The following is another example on how to use the ImplicitOptimizer describing a biological autoregulation process using two coupled implicit equations.

using DataDrivenDiffEq
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using LinearAlgebra
using DataDrivenSparse
using Plots
gr()

@parameters α = 1.0 β = 1.3 γ = 2.0 δ = 0.5
@variables (x(t))[1:2] = [20.0, 12.0]
eqs = [
    D(x[1]) ~ α / (1 + x[2]) - β * x[1],
    D(x[2]) ~ γ / (1 + x[1]) - δ * x[2],
]
@named autoregulation = System(eqs, t)
sys = mtkcompile(autoregulation)

tspan = (0.0, 5.0)
de_problem = ODEProblem(sys, [], tspan)
de_solution = solve(de_problem, Tsit5(), saveat = 0.005);
plot(de_solution)

As always, we start by defining a DataDrivenProblem and a sufficient basis for sparse regression.

dd_prob = DataDrivenProblem(de_solution)

x = sys.x
eqs = [polynomial_basis(x, 4); D.(x); x .* D(x[1]); x .* D(x[2])]

basis = Basis(eqs, x, independent_variable = t, implicits = D.(x))

plot(dd_prob)

Next to varying over different sparsity penalties, we also want to batch our data using DataProcessing.

sampler = DataProcessing(split = 0.8, shuffle = true, batchsize = 30)
res = solve(
    dd_prob, basis, ImplicitOptimizer(STLSQ(1.0e-2:1.0e-2:1.0)),
    options = DataDrivenCommonOptions(data_processing = sampler, digits = 2)
)

"DataDrivenSolution{Float64}" with 2 equations and 10 parameters.
Returncode: Success
Residual sum of squares: 7.788390727504896

And have a look at the resulting plot

plot(res)

As well as the resulting equations of motion

Model ##Basis#435 with 2 equations
States : (x(t))[1] (x(t))[2]
Parameters : 10
Independent variable: t
Equations
Differential(t, 1)((x(t))[1]) = (p₁ + p₂*(x(t))[1] + p₃*(x(t))[1]*(x(t))[2]) / (-p₄ - p₅*(x(t))[2])
Differential(t, 1)((x(t))[2]) = (p₆ + p₇*(x(t))[2] + p₈*(x(t))[1]*(x(t))[2]) / (-p₉ - p₁₀*(x(t))[1])

Have a look at the parameter values

get_parameter_values(system)

10-element Vector{Float64}:
  2.0
 -0.5
 -0.5
 -1.0
 -1.0
  0.76
 -1.0
 -1.0
 -0.76
 -0.76

And the coefficient of determination of the result.

r2(res)

0.9998005648635765

Copy-Pasteable Code

using DataDrivenDiffEq
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using LinearAlgebra
using DataDrivenSparse

@parameters α = 1.0 β = 1.3 γ = 2.0 δ = 0.5
@variables (x(t))[1:2] = [20.0, 12.0]
eqs = [
    D(x[1]) ~ α / (1 + x[2]) - β * x[1],
    D(x[2]) ~ γ / (1 + x[1]) - δ * x[2],
]
@named autoregulation = System(eqs, t)
sys = mtkcompile(autoregulation)

tspan = (0.0, 5.0)
de_problem = ODEProblem(sys, [], tspan)
de_solution = solve(de_problem, Tsit5(), saveat = 0.005);

dd_prob = DataDrivenProblem(de_solution)

x = sys.x
eqs = [polynomial_basis(x, 4); D.(x); x .* D(x[1]); x .* D(x[2])]

basis = Basis(eqs, x, independent_variable = t, implicits = D.(x))

sampler = DataProcessing(split = 0.8, shuffle = true, batchsize = 30)
res = solve(
    dd_prob, basis, ImplicitOptimizer(STLSQ(1.0e-2:1.0e-2:1.0)),
    options = DataDrivenCommonOptions(data_processing = sampler, digits = 2)
)

system = get_basis(res) #hide

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

This page was generated using Literate.jl.