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()

@mtkmodel Autoregulation begin
    @parameters begin
        α = 1.0
        β = 1.3
        γ = 2.0
        δ = 0.5
    end
    @variables begin
        (x(t))[1:2] = [20.0; 12.0]
    end
    @equations begin
        D(x[1]) ~ α / (1 + x[2]) - β * x[1]
        D(x[2]) ~ γ / (1 + x[1]) - δ * x[2]
    end
end

@mtkbuild sys = Autoregulation()
tspan = (0.0, 5.0)
de_problem = ODEProblem{true, SciMLBase.NoSpecialize}(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 the TEXT

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

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

And have a look at the resulting plot

plot(res)

As well as the resulting equations of motion

Model ##Basis#371 with 2 equations
States : (x(t))[1] (x(t))[2]
Parameters : 10
Independent variable: t
Equations
Differential(t)((x(t))[1]) = (p₁ + p₂*(x(t))[1] + p₃*(x(t))[1]*(x(t))[2]) / (-p₄ - p₅*(x(t))[2])
Differential(t)((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}:
  0.76
 -0.99
 -1.0
 -0.76
 -0.76
 -1.0
  0.25
  0.24
  0.49
  0.49

And the coefficient of determination of the result.

r2(res)

0.9997800796994417

Copy-Pasteable Code

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

@mtkmodel Autoregulation begin
    @parameters begin
        α = 1.0
        β = 1.3
        γ = 2.0
        δ = 0.5
    end
    @variables begin
        (x(t))[1:2] = [20.0; 12.0]
    end
    @equations begin
        D(x[1]) ~ α / (1 + x[2]) - β * x[1]
        D(x[2]) ~ γ / (1 + x[1]) - δ * x[2]
    end
end

@mtkbuild sys = Autoregulation()
tspan = (0.0, 5.0)
de_problem = ODEProblem{true, SciMLBase.NoSpecialize}(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(1e-2:1e-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.