ODE Solvers, Process Error, and Difference Equations

ChangeComponent() classes

Model specs contain a set of three lists: * before : instructions to run before the simulation loop begins. * during : instructions to run at each iteration of the simulation loop. * after : instructions to run after the simulation loop ends.

The simplest way to define these lists is using two-sided R formulas. But in the during list one may specify different types of list components, which are of class ChangeComponent. Internally, standard R formulas are converted into an object of class Formula, which is the simplest kind of ChangeComponent. Another important ChangeComponent is the PerCapitaFlow type, which is used for standard flows from one compartment to another.

The list of change components, which is just the during field but ensuring that all elements are valid ChangeComponent objects, can be obtained using the change_components() method in a model spec object.

cc = sir$change_components()
print(cc)
#> [[1]]
#> N ~ S + I + R
#> 
#> [[2]]
#> From: S
#> To: I
#> Per-capita rate expression: beta * I/N
#> Absolute rate symbol: infection
#> [[3]]
#> From: I
#> To: R
#> Per-capita rate expression: (1 - p) * gamma
#> Absolute rate symbol: recovery
#> [[4]]
#> From: I
#> To: D
#> Per-capita rate expression: p * gamma
#> Absolute rate symbol: death

All ChangeComponent objects must have the following methods.

cc[[2]]$change_frame()
#>   state     change
#> 1     S -infection
#> 2     I +infection

cc[[1]]$change_frame()
#> [1] state  change
#> <0 rows> (or 0-length row.names)

cc[[2]]$flow_frame()
#>   size    change       rate         abs_rate
#> 1    S infection beta * I/N S * (beta * I/N)

cc[[1]]$flow_frame()
#> [1] size     change   rate     abs_rate
#> <0 rows> (or 0-length row.names)

ChangeModel() classes

The ChangeModel objects combine the information in lists of ChangeComponent objects, so that they specify a full model. It expands the concept of before, during, and after into a more refined set of steps. All of these methods return lists of two-sided formulas.

• before_loop() : The before list formulas.
• once_start() : Formulas to evaluate at the start of the during list, and which will not be repeated (or modified and repeated) throughout expansions of the during loop. An example of such an expansion is a Runge-Kutta differential equation solver that reuses expressions in an iterative manner. Formulas returned by once_during() are useful for specifying exogeneous changes (e.g., time-varying parameters) that occur once per time-step but that should not change throughout the within-time-step iterations of a differential equation solver.
• before_flows() : Formulas to evaluate before the absolute flow rates are computed.
• update_flows() : Formulas that update the flows, using flow_frame() methods of individual ChangeComponent objects.
• before_state() : Formulas to evaluate before the state vector is updated.
• update_state() : Formulas that update the state vector, using change_frame() method s of individual ChangeComponent objects.
• after_state() : Formulas to evaluate after the state vector is updated.
• once_end() : Formulas to evaluate at the end of the during list, and which will not be repeated throughout (similar to once_start()).
• after_loop() : The after list formulas.

The ChangeModel objects also have flow_frame() and change_frame() methods for combining the outputs of these methods within the individual ChangeComponent objects. An example for an SIR model gives the flow_frame() and change_frame() output as follows.

state   change
-----   ------
S       -infection
I       +infection
I       -recovery
R       +recovery

size  change     rate
----  ------     ----
S     infection  beta * I / N
I     recovery   gamma

The update_flows() method The update_state() method makes use of this flow_frame() to produce the

si = mp_tmb_library("starter_models", "sir", package = "macpan2")
si$change_model
#> Classes 'SimpleChangeModel', 'ChangeModelDefaults', 'ChangeModel', 'Base' <environment: 0x558c2e1dae18> 
#> after_loop : function ()  
#> after_state : function ()  
#> all_user_aware_names : function ()  
#> before_flows : function ()  
#> before_loop : function ()  
#> before_state : function ()  
#> change_classes : function ()  
#> change_frame : function ()  
#> flow_frame : function ()  
#> no_change_components : function ()  
#> once_finish : function ()  
#> once_start : function ()  
#> other_generated_formulas : function ()  
#> update_flows : function ()  
#> update_state : function ()  
#> user_formulas : function ()
si$change_model$flow_frame()
#>   size    change       rate         abs_rate
#> 1    S infection I * beta/N S * (I * beta/N)
#> 2    I  recovery      gamma      I * (gamma)
si$change_model$change_frame()
#>   state     change
#> 1     S -infection
#> 2     I +infection
#> 3     I  -recovery
#> 4     R  +recovery
si$change_model$update_state()
#> [[1]]
#> S ~ -infection
#> <environment: 0x558c2d161e60>
#> 
#> [[2]]
#> I ~ +infection - recovery
#> <environment: 0x558c2d161e60>
#> 
#> [[3]]
#> R ~ +recovery
#> <environment: 0x558c2d161e60>
si$change_model$update_flows()
#> $S
#> $S[[1]]
#> infection ~ I * beta/N
#> <environment: 0x558c2cd52320>
#> 
#> 
#> $I
#> $I[[1]]
#> recovery ~ gamma
#> <environment: 0x558c2cd52320>

spec = mp_tmb_model_spec(
  during = list(
      mp_per_capita_flow("A", "B", "a", "r")
    , mp_per_capita_flow("B", "C", "b", "rr")
  )
)
spec$change_model$update_flows()
#> $A
#> $A[[1]]
#> r ~ a
#> <environment: 0x558c2b5c6770>
#> 
#> 
#> $B
#> $B[[1]]
#> rr ~ b
#> <environment: 0x558c2b5c6770>

UpdateMethod() classes

In order to define a state updater one must define a new UpdateMethod class. These classes are required to update three methods that do not take any arguments: before, during, and after. Each of these methods are required to return a list of two-sided formulas giving the expression list for the three phases of the simulation: before the simulation loop, at every iteration of the simulation loop, and after the simulation loop. Examples of UpdateMethod classes include: EulerUpdateMethod, RK4UpdateMethod, EulerMultinomialUpdateMethod, and HazardUpdateMethod.

The constructors of these UpdateMethod classes often have a field of class ChangeModel, which specifies how to sort the components of

Euler

The simplest approach is to just pretend that the ODEs are difference equations. In this case, at each time-step, each state variable is updated as follows. x_i = x_i − ∑_jx_ir_ij + ∑_jx_jr_ji The SIR example is as follows.

S = S − Sr_SI I = I − Ir_IR + Sr_SI R = R + Ir_IR

Runge Kutta 4

TODO

Euler-Multinomial

The Euler-multinomial state update method assumes that the number of individuals that move from one box to another in a single time-step is a random non-negative integer, coming from a particular multinomial distribution that we now describe.

The probability of staying in the ith box through an entire time-step is assumed to be the following (TODO: relate this to Poisson processes).

p_ii = exp (−∑_jr_ij)

This probability assumes that the r_ij are held constant throughout the entire time-step, although they can change when each time-step begins.

The probability of moving from box i to box j in one time-step is given by the following.

$$ p_{ij} = (1 - p_{ii}) \frac{r_{ij}}{\sum_j r_{ij}} $$

This probability is just the probability of not staying in box i, which is 1 − p_ii, and specifically going to box j, which is assumed to be given by this expression $\frac{r_{ij}}{\sum_j r_{ij}}$.

Let ρ_ij be the random number of individuals that move from box i to box j in one time-step. The expected value of ρ_ij is p_ijx_i. However, these random variables are not independent events, because the total number of individuals, ∑_ix_i, has to remain constant through a single time-step (at least in the models that we are currently considering).

To account for this non-independence, we collect the ρ_ij associated with a from compartment i into a vector ρ_i. We collect similar vector of probabilities, p_i. Each ρ_i is a random draw from a multinomial distribution with x_i trials and probability vector, p_i.

Once these random draws have been made, the state variables can be updated at each time-step as follows.

x_i = x_i − ∑_jρ_ij + ∑_jρ_ji

Note that we do not actually need to generate values for the diagonal elements, ρ_ii, because they cancel out in this update equation. We also can ignore any ρ_ij such that r_ij = 0.

Under the SIR example we have a particularly simple Euler-binomial distribution because there are no branching flows – when individuals leave S they can only go to I and when they leave I they can only go to R. These two flows are given by the following distributions.

ρ_SI ∼ Binomial(S, p_SI)

ρ_IR ∼ Binomial(I, p_IR)

In models with branching flows we would have multinomial distributions. But in this model the state update is given by the following equations.

S = S − ρ_SI I = I − ρ_IR + ρ_SI R = R + ρ_IR

Hazard

The Hazard update-method uses the expected values associated with the Euler-multinomial described above. In particular, the state variables are updated as follows.

x_i = x_i − ∑_jx_ip_ij + ∑_jx_jp_ji

The SIR model would be as follows.

S = S − Sp_SI I = I − Ip_IR + Sp_SI R = R + Ip_IR

More explicitly for the I compartment this would be. I(t + 1) = I(t) − I(t)(1 − exp (−γ)) + S(t)(1 − exp (−βI(t)))

Linearizing at Each Time-Step

Because the r_ij could depend on any state variable, the ODE above is generally non-linear. However, we could linearize the model at every time-step and explicitly compute the matrix exponential to find the approximate state update.

Hazard in models including more than unbalanced per-capita flows

The perfectly balanced per-capita flows approach used above does not always work. For example, virus shedding to a wastewater compartment does not involve infected individuals flowing into a wastewater compartment, because people do not become wastewater obviously. In such models we need to think more clearly about how to use the Hazard approximation.

For such cases we can modify the differential equation to include an arbitrary number of absolute inflows and outflows.

$$ \frac{dx_i}{dt} = \underbrace{\sum_j x_j r_{ji}}_{\text{inflow}} - \underbrace{\sum_j x_i r_{ij}}_{\text{outflow}} + \underbrace{\sum_k r^+_{ik}}_{\text{absolute inflow}} - \underbrace{\sum_l r^-_{il}}_{\text{absolute outflow}} $$

Where r_ik⁺ is the absolute rate at which x_i increases due to mechanism k and r_il⁻ is the absolute rate at which x_i decreases due to mechanism l.