Spline Example: A Modulating Input Function
In this example, it is demonstrated how a continuous input function that modulates the right-hand side of an ODE in a time-dependent manner can be estimated without direct observations. The model used for demonstration purpose is an SEIR model. These models are frequently applied for describing the dynamics of infection diseases. In this context, the time-dependent modulation originates from political containment actions that reduces the infection rate.
In particular, instead of a commont rate or reaction
Susceptibles -> Exposed CUSTOM "b * Infectious * Susceptibles" "description"
it might be desired to have time dependent b, e.g.
Susceptibles -> Exposed CUSTOM "b * b_time_dependence * Infectious * Susceptibles" "description"
where b_time_dependence denotes a cumstom function of time.
The aim of this example is not having a realistic model and a common parametrization from literature. Instead, the model is rather simple for demonstrating capabilities of D2D in input estimation.
The implementation is available in folder Examples/ToyModels/SplineExample_InputModulation
In the example, the infection rate is known to be time dependent due to containmenta actions but the functional form is unknown. Moreover, the infection rate is not observed directly. Only the infected people termed "Infectious" are observed.
Smoothing splines are frequently applied for estimating smooth functional relationships. Splines are polynomials that are define in intervals and connected a the "knots" via boundary conditions that ensure a continuous up to the 2nd derivative.
The shape of a smoothing spline depends on the data and the degrees of freedoms, i.e. the number of knots. In our implementation, one knot corresponds to one parameter. A parametrization is chosen that the parameter depends on the value of the spline at each knot.
In the context of ODE-modelling, one typical challenge of using splines is that negative values should be prevented. This can be done by using the the spline as exponent, i.e. the input is given by 10^spline(t). In our context, the restrictions are even stronger because the spline should be within the range (0,1]. Moreover, it was intended for prediction purpose to switch-off the impact of the smoothing spline after the last observation.
In our implemenation, the smoothing spline s(t) with 10 knots was linked with a sigmoidal function
f(s) = 0.5*(1+(1-exp(-s))/(1+exp(-s)))*(1-s)
This function f is zero for small s and one for large s.
In the model definition file, the implementation reads
INPUTS
InputSpline C 1 1 "spline10(t, 0.0, 5, 10.0, sp2, 20.0, sp3, 30.0, sp4, 40.0, sp5, 50.0, sp6, 60.0, sp7, 80.0, sp8, 150.0, sp9, 400.0, 5, 0, 0.0)"
...
DERIVED
b_time_dependence C "N" "N" "0.5*(1+(1-exp(-InputSpline))/(1+exp(-InputSpline)))"
Note, that the spline has ten knots, at t=0, 10, 20, 30, 40, 50, 60, 80, 150, 400 but only 8 free parameters. The first and the last spline parameters are set to 5 because f(5) is close to one. This was a further intended constaint of the model.
In order to keep the input constant after the last data point, a second spline-input is defined that shares the parameters sp2, ..., sp9. Instead of the model time t, the second spline depends on the argument tLast. This means, that the second spline has the same response as the first for t=tLast. In order to switch from the first to the second spline at t=tLast, both splines are connected with an a step input in the following way:
INPUTS
InputSpline C 1 1 "spline10(t, 0.0, 5, 10.0, sp2, 20.0, sp3, 30.0, sp4, 40.0, sp5, 50.0, sp6, 60.0, sp7, 80.0, sp8, 150.0, sp9, 400.0, 5, 0, 0.0)"
InputSplineEnd C 1 1 "spline10(tLast, 0.0, 5, 10.0, sp2, 20.0, sp3, 30.0, sp4, 40.0, sp5, 50.0, sp6, 60.0, sp7, 80.0, sp8, 150.0, sp9, 400.0, 5, 0, 0.0)"
fixSpline C 1 1 "step1(t, 0, tLast, 1)"
...
DERIVED
b_time_dependence C "N" "N" "0.5*(1+(1-exp(-InputSpline))/(1+exp(-InputSpline)))*(1-fixSpline) + fixSpline*0.5*(1+(1-exp(-InputSplineEnd))/(1+exp(-InputSplineEnd)))"