Create SPR_fitting.md

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))
+```

src/reaction_network.jl

@@ -15,7 +15,7 @@
 @reaction_network begin
     p, 0 --> X
     d, X --> 0
-end p d
+end 
 ```
 
 Parameters and species can be explicitly indicated using the @parameters and @species
@@ -123,7 +123,7 @@
 end
 
 # Declares the keys used for various options.
-const option_keys = (:species, :parameters, :variables, :ivs)
+const option_keys = (:species, :parameters, :variables, :ivs, :discrete_events, :continuous_events)
 
 ### The @species macro, basically a copy of the @variables macro. ###
 macro species(ex...)
@@ -156,6 +156,14 @@
     esc(vars)
 end
 
+macro discrete_events(ex...)
+    ex
+end
+
+macro continuous_events(ex...)
+    ex
+end
+
 ### The main macro, takes reaction network notation and returns a ReactionSystem. ###
 """
     @reaction_network
@@ -354,15 +362,20 @@
     reaction_lines = Expr[x for x in ex.args if x.head == :tuple]
     option_lines = Expr[x for x in ex.args if x.head == :macrocall]
 
+    @show option_lines
+
     # Get macro options.
     options = Dict(map(arg -> Symbol(String(arg.args[1])[2:end]) => arg,
                        option_lines))
 
+    @show options
     # Parses reactions, species, and parameters.
     reactions = get_reactions(reaction_lines)
     species_declared = extract_syms(options, :species)
     parameters_declared = extract_syms(options, :parameters)
     variables = extract_syms(options, :variables)
+    discrete_e = extract_discrete_events(options)
+    continuous_e = extract_continuous_events(options)
 
     # handle independent variables
     if haskey(options, :ivs)
@@ -413,7 +426,9 @@
 
         Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, union($spssym, $varssym),
                                               $pssym; name = $name,
-                                              spatial_ivs = $sivs)
+                                              spatial_ivs = $sivs,
+                                              discrete_events = $discrete_e,
+                                              continuous_events = $continuous_e)
     end
 end
 
@@ -474,6 +489,51 @@
     return syms
 end
 
+# When the user have used the @continuous_events option, extract continuous events
+function extract_continuous_events(opts)
+    events = Vector{Equation}[]
+    if haskey(opts, :continuous_event)
+        push!(events, opts[:continuous_events])
+    end
+    if length(events)==0
+        return nothing
+    else
+        return events
+    end
+end
+
+# When the user have used the @discrete_events option, extract discrete events
+#function extract_discrete_events(opts)
+    #events = Expr[]
+    #if haskey(opts, :discrete_events)
+        #push!(events, opts[:discrete_events])
+    #end
+    #if length(events)==0
+        #return nothing
+    #else
+        #return events
+    #end
+#end
+
+#this is interpolating rather than quoting. why does it need to be quoting? 
+function extract_discrete_events(opts)
+    events = []
+    if haskey(opts, :discrete_events)
+        ex = quote 
+            $(opts[:discrete_events])
+        end
+        ex1 = MacroTools.striplines(ex)
+        ex2 = ex1.args[1].args[end].args
+        push!(events, ex2)
+    end
+    if length(events)==0
+        return nothing
+    else
+        return events
+    end
+end
+
+#interp
 # Function looping through all reactions, to find undeclared symbols (species or
 # parameters), and assign them to the right category.
 function extract_species_and_parameters!(reactions, excluded_syms)