|
| 1 | +# [Inputs and time-dependent (or functional) parameters](@id time_dependent_parameters) |
| 2 | +Catalyst supports the usage of "functional parameters". In practice, these are parameters that are given by (typically) time-dependent functions (they can also depend on e.g. species values, as discussed [here](@ref functional_parameters_sir)). They are a way to inject custom functions into models. Functional parameters can be used when rates depend on real data, or to represent complicated functions (which use e.g. `for` loops or random number generation). Here, the function's values are declared as a data interpolation (which interpolates discrete samples to a continuous function). This is then used as the functional parameter's value in the simulation. This tutorial first shows how to create time-dependent functional parameters, and then gives an example where the functional parameter depends on a species value. |
| 3 | + |
| 4 | +An alternative approach for representing complicated functions is by [using `@register_symbolic`](@ref dsl_description_nonconstant_rates_function_registration). |
| 5 | + |
| 6 | +## [Basic example](@id functional_parameters_basic_example) |
| 7 | +Let us first consider an easy, quick-start example (the next section will discuss what is going on in more detail). We will consider a simple [birth-death model](@ref basic_CRN_library_bd), but where the birth rate is determined by an input parameter (for which the value depends on time). First, we [define the input parameter programmatically](@ref programmatic_CRN_construction), and its values across all time points using the [DataInterpolations.jl](https://github.com/SciML/DataInterpolations.jl) package. In this example we will use the input function $pIn(t) = (2 + t)/(1 + t)$. Finally, we plot the input function, demonstrating how while it is defined at discrete points, DataInterpolations.jl generalises this to a continuous function. |
| 8 | +```@example functional_parameters_basic_example |
| 9 | +using Catalyst, DataInterpolations, Plots |
| 10 | +t = default_t() |
| 11 | +tend = 10.0 |
| 12 | +ts = collect(0.0:0.01:tend) |
| 13 | +spline = LinearInterpolation((2 .+ ts) ./ (1 .+ ts), ts) |
| 14 | +@parameters (pIn::typeof(spline))(..) |
| 15 | +plot(spline) |
| 16 | +``` |
| 17 | +Next, we create our model, [interpolating](@ref dsl_advanced_options_symbolics_and_DSL_interpolation) the input parameter into it (making it a function of `t`). |
| 18 | +```@example functional_parameters_basic_example |
| 19 | +input = pIn(default_t()) |
| 20 | +bd_model = @reaction_network begin |
| 21 | + $input, 0 --> X |
| 22 | + d, X --> 0 |
| 23 | +end |
| 24 | +``` |
| 25 | +Finally, we can simulate our model as normal (but where we set the value of the `pIn` parameter to our interpolated data). |
| 26 | +```@example functional_parameters_basic_example |
| 27 | +using OrdinaryDiffEqDefault |
| 28 | +u0 = [:X => 0.5] |
| 29 | +ps = [:d => 2.0, :pIn => spline] |
| 30 | +oprob = ODEProblem(bd_model, u0, tend, ps) |
| 31 | +sol = solve(oprob) |
| 32 | +plot(sol) |
| 33 | +``` |
| 34 | +!!! note |
| 35 | + For this simple example, $(2 + t)/(1 + t)$ could have been used directly as a reaction rate (or written as a normal function), technically making the functional parameter approach unnecessary. However, here we used this function as a simple example of how discrete data can be made continuous using DataInterpolations, and then have its values inserted using a (functional) parameter. |
| 36 | + |
| 37 | +## [Inserting a customised, time-dependent, input](@id functional_parameters_circ_rhythm) |
| 38 | +Let us now go through everything again, but providing some more details. Let us first consider the input parameter. We have previously described how a [time-dependent rate can model a circadian rhythm](@ref dsl_description_nonconstant_rates_time). For real applications, due to e.g. clouds, sunlight is not a perfect sine wave. Here, a common solution is to take real sunlight data from some location and use in the model. Here, we will create synthetic (noisy) data as our light input: |
| 39 | +```@example functional_parameters_circ_rhythm |
| 40 | +using Plots |
| 41 | +tend = 120.0 |
| 42 | +ts = collect(0.0:1.0:tend) |
| 43 | +light = sin.(ts/6) .+ 1 |
| 44 | +light = [max(0.0, l - rand()) for l in light] |
| 45 | +plot(ts, light; seriestype = :scatter, label = "Experienced light") |
| 46 | +``` |
| 47 | +Now this input is only actually defined at the sample points, making it incompatible with a continuous ODE simulation. To enable this, we will use the DataInterpolations package to create an interpolated version of this data, which forms the actual input: |
| 48 | +```@example functional_parameters_circ_rhythm |
| 49 | +using DataInterpolations |
| 50 | +interpolated_light = LinearInterpolation(light, ts) |
| 51 | +plot(interpolated_light) |
| 52 | +``` |
| 53 | +We are now ready to declare our model. We will consider a protein with an active and an inactive form ($Pₐ$ and $Pᵢ$) where the activation is driven by the presence of sunlight. In this example we we create our model using the [programmatic approach](@ref programmatic_CRN_construction). Do note the special syntax we use to declare our input parameter, where we both designate it as a generic function and its type as the type of our interpolated input. Also note that, within the model, we mark the input parameter (`light_in`) as a function of `t`. |
| 54 | +```@example functional_parameters_circ_rhythm |
| 55 | +using Catalyst |
| 56 | +t = default_t() |
| 57 | +in_type = typeof(interpolated_light) |
| 58 | +@parameters kA kD (light_in::in_type)(..) |
| 59 | +@species Pₐ(t) Pᵢ(t) |
| 60 | +rxs = [ |
| 61 | + Reaction(kA*light_in(t), [Pᵢ], [Pₐ]), |
| 62 | + Reaction(kD, [Pₐ], [Pᵢ]) |
| 63 | +] |
| 64 | +@named rs = ReactionSystem(rxs, t) |
| 65 | +rs = complete(rs) |
| 66 | +``` |
| 67 | +Now we can simulate our model. Here, we use the interpolated data as the input parameter's value. |
| 68 | +```@example functional_parameters_circ_rhythm |
| 69 | +using OrdinaryDiffEqDefault |
| 70 | +u0 = [Pᵢ => 1.0, Pₐ => 0.0] |
| 71 | +ps = [kA => 1.5, kD => 1.0, light_in => interpolated_light] |
| 72 | +oprob = ODEProblem(rs, u0, tend, ps) |
| 73 | +sol = solve(oprob) |
| 74 | +plot(sol) |
| 75 | +``` |
| 76 | + |
| 77 | +### [Interpolating the input into the DSL](@id functional_parameters_circ_rhythm_dsl) |
| 78 | +It is possible to use time-dependent inputs when creating models [through the DSL](@ref dsl_description) as well. However, it can still be convenient to declare the input parameter programmatically as above. Using it, we form an expression of it as a function of time, and then [interpolate](@ref dsl_advanced_options_symbolics_and_DSL_interpolation) it into our DSL-declaration: |
| 79 | +```@example functional_parameters_circ_rhythm |
| 80 | +input = light_in(t) |
| 81 | +rs_dsl = @reaction_network rs begin |
| 82 | + (kA*$input, kD), Pᵢ <--> Pₐ |
| 83 | +end |
| 84 | +``` |
| 85 | +We can confirm that this model is identical to our programmatic one (and should we wish to, we can simulate it using identical syntax). |
| 86 | +```@example functional_parameters_circ_rhythm |
| 87 | +rs == rs_dsl |
| 88 | +``` |
| 89 | + |
| 90 | +## [Non-time functional parameters](@id functional_parameters_sir) |
| 91 | +Previously we have demonstrated functional parameters that are functions of time. However, functional parameters can be functions of any variable (however, currently, more than one argument is not supported). Here we will demonstrate this using a [SIR model](@ref basic_CRN_library_sir), but instead of having the infection rate scale linearly with the number of infected individuals, we instead assume we have measured data of the infection rate (as dependent on the number of infected individuals) and wish to use this instead. Normally we use the following infection reaction in the SIR model: |
| 92 | +```julia |
| 93 | +@reaction k1, S + I --> 2I |
| 94 | +``` |
| 95 | +For ODE models, this would give the same equations as |
| 96 | +```julia |
| 97 | +@reaction k1*I, S --> I |
| 98 | +``` |
| 99 | +Due to performance reasons (especially for jump simulations) the former approach is *strongly* encouraged. Here, however, we will assume that we have measured real data of how the number of infected individuals affects the infection rate, and wish to use this in our model, i.e. something like this: |
| 100 | +```julia |
| 101 | +@reaction k1*i_rate(I), S --> I |
| 102 | +``` |
| 103 | +This is a case where we can use a functional parameter, whose value depends on the value of $I$. |
| 104 | + |
| 105 | +We start by declaring the functional parameter that describes how the infection rate depends on the number of infected individuals. We also plot the measured infection rate, and compare it to the theoretical rate usually used in the SIR model. |
| 106 | +```@example functional_parameters_sir |
| 107 | +using DataInterpolations, Plots |
| 108 | +I_grid = collect(0.0:5.0:100.0) |
| 109 | +I_measured = 300.0 *(0.8*rand(length(I_grid)) .+ 0.6) .* I_grid ./ (300 .+ I_grid) |
| 110 | +I_rate = LinearInterpolation(I_measured, I_grid) |
| 111 | +plot(I_rate; label = "Measured infection rate") |
| 112 | +plot!(I_grid, I_grid; label = "Normal SIR infection rate") |
| 113 | +``` |
| 114 | +Next, we create our model (using the DSL approach). As `I_rate` will be a function of $I$ we will need to declare this species first as well. |
| 115 | +```@example functional_parameters_sir |
| 116 | +using Catalyst |
| 117 | +@parameters (inf_rate::typeof(I_rate))(..) |
| 118 | +@species I(default_t()) |
| 119 | +inf_rate_in = inf_rate(I) |
| 120 | +sir = @reaction_network rs begin |
| 121 | + k1*$inf_rate_in, S --> I |
| 122 | + k2, I --> R |
| 123 | +end |
| 124 | +nothing # hide |
| 125 | +``` |
| 126 | +Finally, we can simulate our model. |
| 127 | +```@example functional_parameters_sir |
| 128 | +using OrdinaryDiffEqDefault |
| 129 | +u0 = [:S => 99.0, :I => 1.0, :R => 0.0] |
| 130 | +ps = [:k1 => 0.002, :k2 => 0.01, :inf_rate => I_rate] |
| 131 | +oprob = ODEProblem(sir, u0, 250.0, ps) |
| 132 | +sol = solve(oprob) |
| 133 | +plot(sol) |
| 134 | +``` |
| 135 | +!!! note |
| 136 | + When declaring a functional parameter of time, it is easy to set its grid values (i.e. ensure they range from the first to the final time point). For Functional parameters that depend on species concentrations it is trickier, and one must make sure that any potential input-species values that can occur during the simulation are represented in the interpolation. |
| 137 | + |
| 138 | +### [Using different data interpolation approaches](@id functional_parameters_interpolation_algs) |
| 139 | +Up until now we have used [linear interpolation](https://en.wikipedia.org/wiki/Linear_interpolation) of our data. However, DataInterpolations.jl [supports other interpolation methods](https://docs.sciml.ai/DataInterpolations/stable/methods/). To demonstrate these we here generate a data set, and then show the linear, cubic, and constant interpolations: |
| 140 | +```@example functional_parameters_interpolation_algs |
| 141 | +using DataInterpolations, Plots |
| 142 | +xs = collect(0.0:1.0:10.0) |
| 143 | +ys = xs ./ (5*rand(length(xs)) .+ xs) |
| 144 | +spline_linear = LinearInterpolation(ys, xs) |
| 145 | +spline_cubuc = CubicSpline(ys, xs) |
| 146 | +spline_const = ConstantInterpolation(ys, xs) |
| 147 | +plot(spline_linear) |
| 148 | +plot!(spline_cubuc) |
| 149 | +plot!(spline_const) |
| 150 | +``` |
| 151 | +Finally, DataInterpolations.jl also allows various [extrapolation methods](https://docs.sciml.ai/DataInterpolations/stable/extrapolation_methods/). |
0 commit comments