From 1945c5faad6e20e6cf5f799a5f719f2619cf9cd3 Mon Sep 17 00:00:00 2001 From: sophie strey <106631595+tstrey@users.noreply.github.com> Date: Wed, 25 Oct 2023 12:19:29 -0400 Subject: [PATCH 1/7] Create SPR_fitting.md --- docs/src/example_networks/SPR_fitting.md | 105 +++++++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 docs/src/example_networks/SPR_fitting.md diff --git a/docs/src/example_networks/SPR_fitting.md b/docs/src/example_networks/SPR_fitting.md new file mode 100644 index 0000000000..413a4cfbc8 --- /dev/null +++ b/docs/src/example_networks/SPR_fitting.md @@ -0,0 +1,105 @@ +# [Surface Plasmon Resonance Model Example](@id surface_plasmon_resonance_model_Example) +This tutorial shows how to programmatically construct a ['ReactionSystem'](@ref) corresponding to the chemistry underlying the [Surface Plasmon Resonance Model](https://en.wikipedia.org/wiki/Surface_plasmon_resonance) using [ModelingToolkit](http://docs.sciml.ai/ModelingToolkit/stable/)/[Catalyst](http://docs.sciml.ai/Catalyst/stable/). You can find a simpler constructruction of this model [here](https://docs.sciml.ai/Catalyst/stable/catalyst_applications/parameter_estimation/) in the example of parameter estimation. + +Using the symbolic interface, we will create our states/species, define our reactions, and add a discrete event. This model is simulated and used to generate sample data points with added noise and then used to fit parameters in iterations. This event will correspond to the time at which the antibody stop binding to antigen and switch to dissociating. + +We begin by importing some necessary packages. +```julia +using ModelingToolkit, Catalyst, DifferentialEquations +using Plots, Optimization +using OptimizatizationOptimJL, Optim +``` +In this example, the concentration of antigen,`\beta` is varied to determine the constant of proportionality, `\alpha`, `k_{on}`(association rate constant), and `k_{off}` (dissociation rate constant) which characterized the binding interaction between the antigen and the immobilized receptor molecules on the sensor slide. We start by defining our reaction equations, parameters, variables, event, and states. +```julia +@variables t +@parameters k_on k_off α +@species A(t) B(t) + +rxs = [(@reaction α*k_on, A --> B), (@reaction k_off, B --> A)] + +switch_time = 2.0 +discrete_events = (t == switch_time) => [k_on ~ 0.0] + +tspan = (0.0, 4.0) +alpha_list = [0.1,0.2,0.3,0.4] #list of concentrations +``` +Iterating over values of `\alpha`, now we create a list of ODE solutions and set the initial conditions of our states and parameters. +```julia +results_list = [] + +u0 = [:A => 10.0, :B => 0.0] +p_real = [k_on => 100.0, k_off => 10.0, α => 1.0] +@named osys = ReactionSystem(rxs, t, [A, B], [k_on, k_off, α]; discrete_events) +oprob = ODEProblem(osys, u0, tspan, p_real) +sample_times = range(tspan[1]; stop = tspan[2], length = 1001) + +for alpha in alpha_list + p_real = [k_on => 100.0, k_off => 10.0, α => alpha] + oprobr = remake(oprob, p=p_real) + sol_real = solve(oprobr, Tsit5(); tstops = sample_times) + + push!(results_list, sol_real(sample_times)) +end +``` +```@example ceq3 +default(; lw = 3, framestyle = :box, size = (800, 400)) +p = plot() +plot(p, sample_times, results_list[1][2,:]) +plot!(p, sample_times, results_list[2][2,:]) +plot!(p, sample_times, results_list[3][2,:]) +plot!(p, sample_times, results_list[4][2,:]) + +sample_vals = [] +for result in results_list + sample_val = Array(result) + sample_val .*= (1 .+ .1 * rand(Float64, size(sample_val)) .- .01) + push!(sample_vals, sample_val) +end + +for val in sample_vals + scatter!(p, sample_times, val[2,:]; color = [:blue :red], legend = nothing) +end +plot(p) +``` +Next, we create a function to fit the parameters. Here we are using `NelderMead()`. In order to achieve a better fit, we are incorporating all of our solutions into the loss function. We will fit seperately the association and dissociation signals so for the first estimate, `tend < switch_time`. +```julia +function optimise_p(pinit, tend) + function loss(p, _) + newtimes = filter(<=(tend), sample_times) + solutions = [] + for alpha in alpha_list + newprob = remake(oprob; tspan = (0.0, tend), p = [k_on => p[1],k_off => p[2],α => alpha]) + sol = Array(solve(newprob, Tsit5(); saveat = newtimes, tstops = switch_time)) + push!(solutions,sol[2,:]) + end + loss = 0 + for (idx, solution) in enumerate(solutions) + loss += sum(abs2, p[3]*solution .- sample_vals[idx][2, 1:size(sol,2)]) + end + + return loss + end + + optf = OptimizationFunction(loss) + + optprob = OptimizationProblem(optf, pinit) + sol = solve(optprob, Optim.NelderMead()) + + return sol.u +end + +p_estimate = optimise_p([100.0, 10.0, 1.0], 1.5) +``` +Finally, we remake the solution using the estimate, fit the entire time span for some `alpha`, and plot the solution. +```julia +alpha = 0.2 +newprob = remake(oprob; tspan = (0.0, 4.0), p = [k_on => p_estimate[1], k_off => p_estimate[2], α => alpha]) +newsol = solve(newprob, Tsit5(); tstops = switch_time) + +#plotting simulated sample date with new solution +default(; lw = 3, framestyle = :box, size = (800, 400)) + +scatter!(sample_times, sample_vals; color = [:darkblue :darkred], legend = nothing) + +plot!(newsol[2,:]; legend = nothing, color = [:blue :red], linestyle = :dash, tspan= (0.0, 4.0)) +``` From db26955537ed90a5f56818729d7ff593f2df13cd Mon Sep 17 00:00:00 2001 From: sophie strey <106631595+tstrey@users.noreply.github.com> Date: Thu, 21 Dec 2023 18:50:28 -0500 Subject: [PATCH 2/7] changing optimization package combination --- docs/src/example_networks/SPR_fitting.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/example_networks/SPR_fitting.md b/docs/src/example_networks/SPR_fitting.md index 413a4cfbc8..2fd16f385c 100644 --- a/docs/src/example_networks/SPR_fitting.md +++ b/docs/src/example_networks/SPR_fitting.md @@ -6,8 +6,8 @@ Using the symbolic interface, we will create our states/species, define our reac We begin by importing some necessary packages. ```julia using ModelingToolkit, Catalyst, DifferentialEquations -using Plots, Optimization -using OptimizatizationOptimJL, Optim +using Optimization, OptimizatizationOptimJL +using Plots ``` In this example, the concentration of antigen,`\beta` is varied to determine the constant of proportionality, `\alpha`, `k_{on}`(association rate constant), and `k_{off}` (dissociation rate constant) which characterized the binding interaction between the antigen and the immobilized receptor molecules on the sensor slide. We start by defining our reaction equations, parameters, variables, event, and states. ```julia From b3c6f93031e35895a5004ca8e84af7769dd94770 Mon Sep 17 00:00:00 2001 From: "Helmut H. Strey" Date: Tue, 5 Nov 2024 10:28:23 -0500 Subject: [PATCH 3/7] remove differential_equations import and change model into DSL --- docs/pages.jl | 3 ++- docs/src/example_networks/SPR_fitting.md | 21 +++++++++++++-------- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/docs/pages.jl b/docs/pages.jl index ca2241c5b7..fe538e7af5 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -14,6 +14,7 @@ pages = Any["Home" => "index.md", "catalyst_applications/parameter_estimation.md"], "Example Networks" => Any["example_networks/basic_CRN_examples.md", "example_networks/hodgkin_huxley_equation.md", - "example_networks/smoluchowski_coagulation_equation.md"], + "example_networks/smoluchowski_coagulation_equation.md", + "example_networks/SPR_fitting.md"], "FAQs" => "faqs.md", "API" => "api/catalyst_api.md"] diff --git a/docs/src/example_networks/SPR_fitting.md b/docs/src/example_networks/SPR_fitting.md index 2fd16f385c..83f8e23b0c 100644 --- a/docs/src/example_networks/SPR_fitting.md +++ b/docs/src/example_networks/SPR_fitting.md @@ -5,23 +5,29 @@ Using the symbolic interface, we will create our states/species, define our reac We begin by importing some necessary packages. ```julia -using ModelingToolkit, Catalyst, DifferentialEquations +using ModelingToolkit, Catalyst using Optimization, OptimizatizationOptimJL using Plots ``` In this example, the concentration of antigen,`\beta` is varied to determine the constant of proportionality, `\alpha`, `k_{on}`(association rate constant), and `k_{off}` (dissociation rate constant) which characterized the binding interaction between the antigen and the immobilized receptor molecules on the sensor slide. We start by defining our reaction equations, parameters, variables, event, and states. ```julia -@variables t -@parameters k_on k_off α -@species A(t) B(t) +osys = @reaction_network begin + @variables t + @parameters k_on=100.0 k_off α + @species A(t)B(t) + + @discrete_events begin + t == switch_time => [k_on ~ 0.0] + end -rxs = [(@reaction α*k_on, A --> B), (@reaction k_off, B --> A)] + α*k_on, A --> B + k_off, B --> A +end switch_time = 2.0 -discrete_events = (t == switch_time) => [k_on ~ 0.0] tspan = (0.0, 4.0) -alpha_list = [0.1,0.2,0.3,0.4] #list of concentrations +alpha_list = [0.1, 0.2, 0.3, 0.4] #list of concentrations ``` Iterating over values of `\alpha`, now we create a list of ODE solutions and set the initial conditions of our states and parameters. ```julia @@ -29,7 +35,6 @@ results_list = [] u0 = [:A => 10.0, :B => 0.0] p_real = [k_on => 100.0, k_off => 10.0, α => 1.0] -@named osys = ReactionSystem(rxs, t, [A, B], [k_on, k_off, α]; discrete_events) oprob = ODEProblem(osys, u0, tspan, p_real) sample_times = range(tspan[1]; stop = tspan[2], length = 1001) From 8223bd2646072d201f2361120c35b8e937b503eb Mon Sep 17 00:00:00 2001 From: sophie strey <106631595+tstrey@users.noreply.github.com> Date: Thu, 7 Nov 2024 09:45:12 -0500 Subject: [PATCH 4/7] Update pages.jl --- docs/pages.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/pages.jl b/docs/pages.jl index 62b773035c..0bf05ef30d 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -23,7 +23,7 @@ pages = Any[ "model_creation/examples/programmatic_generative_linear_pathway.md", "model_creation/examples/hodgkin_huxley_equation.md", "model_creation/examples/smoluchowski_coagulation_equation.md", - "model_creation/examples/SPR_fitting.md" + "model_creation/examples/SPR_fitting.md"] ], "Model simulation and visualization" => Any[ "model_simulation/simulation_introduction.md", From 04a9b7087d83526cb297b468305bf1d11ef37892 Mon Sep 17 00:00:00 2001 From: sophie strey <106631595+tstrey@users.noreply.github.com> Date: Thu, 7 Nov 2024 10:01:08 -0500 Subject: [PATCH 5/7] Rename docs/src/example_networks/SPR_fitting.md to docs/src/model_creation/examples/SPR_fitting.md --- .../{example_networks => model_creation/examples}/SPR_fitting.md | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename docs/src/{example_networks => model_creation/examples}/SPR_fitting.md (100%) diff --git a/docs/src/example_networks/SPR_fitting.md b/docs/src/model_creation/examples/SPR_fitting.md similarity index 100% rename from docs/src/example_networks/SPR_fitting.md rename to docs/src/model_creation/examples/SPR_fitting.md From ddac3bfb151b41032b655058f7940dc53af9d3d1 Mon Sep 17 00:00:00 2001 From: sophie strey <106631595+tstrey@users.noreply.github.com> Date: Thu, 7 Nov 2024 10:12:36 -0500 Subject: [PATCH 6/7] Rename SPR_fitting.md to SPR_fitting.md --- .../{model_creation => inverse_problems}/examples/SPR_fitting.md | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename docs/src/{model_creation => inverse_problems}/examples/SPR_fitting.md (100%) diff --git a/docs/src/model_creation/examples/SPR_fitting.md b/docs/src/inverse_problems/examples/SPR_fitting.md similarity index 100% rename from docs/src/model_creation/examples/SPR_fitting.md rename to docs/src/inverse_problems/examples/SPR_fitting.md From 9dd92160dd9276e9f8e62f09accfe9e322ed73b6 Mon Sep 17 00:00:00 2001 From: sophie strey <106631595+tstrey@users.noreply.github.com> Date: Thu, 7 Nov 2024 10:13:56 -0500 Subject: [PATCH 7/7] Update pages.jl --- docs/pages.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/docs/pages.jl b/docs/pages.jl index 0bf05ef30d..c03ab89a0a 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -22,8 +22,7 @@ pages = Any[ "model_creation/examples/basic_CRN_library.md", "model_creation/examples/programmatic_generative_linear_pathway.md", "model_creation/examples/hodgkin_huxley_equation.md", - "model_creation/examples/smoluchowski_coagulation_equation.md", - "model_creation/examples/SPR_fitting.md"] + "model_creation/examples/smoluchowski_coagulation_equation.md"] ], "Model simulation and visualization" => Any[ "model_simulation/simulation_introduction.md", @@ -51,8 +50,8 @@ pages = Any[ # "inverse_problems/structural_identifiability.md", "inverse_problems/global_sensitivity_analysis.md", "Examples" => Any[ - "inverse_problems/examples/ode_fitting_oscillation.md" - ] + "inverse_problems/examples/ode_fitting_oscillation.md", + "inverse_problems/examples/SPR_fitting.md"] ], "Spatial modelling" => Any[ "spatial_modelling/lattice_reaction_systems.md",