Fitting a Simple SIR model with parameters that depend on time

Modeling
, ,
1

Hello!

I am trying to fit parameters in ODE-based model using Stan.
The model is from Mathematical modeling of interaction between innate and adaptive immune responses in COVID‐19 and implications for viral pathogenesis, specifically:

image
image837×269 11.7 KB

where the time-dependent parameters are
image

The initial values for T, I, and V are given, and all other parameters are given.

I’ve been able to replicate this successfully in R, however I’m new to Stan and am unsure how to begin inputting this model, especially the time-dependent parametersm in Stan so that I can investigate this study’s choice in parameters givenexperimental data.

Any and all advice are welcome! Thank you!

1 Like
2

Hi. Welcome to the forums. I was going to point you to @charlesm93 (and others)'s excellent case study, but realised people like @wds15, @yizhang and @bbbales2 may know of much better resources than what I recommended. They may also be able to advise on how to implement the piece-wise functions for time-dependent parameters and whether those play nicely with the ODE solvers in Stan.

1 Like
3

Thank you for the response! Yes, I anticipate it will be tricky to get those time-dependent parameters to play nice with the model based on the little I know about Stan.
I’ll be sure to check out the case study you sent and see if that inspires any breakthroughs. If nothing comes of it, I’ll dm the experts you mentioned.

Thanks again!

4

I think it’s best to post publicly so others can learn too. Good luck!

1 Like
5

Hello again!
@wds15, @yizhang and @bbbales2, per @maxbiostat 's suggestion, I am seeking your advice for my question. I’ve given it my best try, and what follows is as far as I’ve gotten on this problem:

functions {
  
  // This largely follows the deSolve package, but also includes the x_r and x_i variables.
  // These variables are used in the background.
  
  real[] deltE(real time,
              real[] peakE,
              real[] wrE,
              real[] wfE) {
                for (i in 1:size(t)){
                  return t<peakE ? exp((-(t-peakE)^2)/wrE) : exp((-(t-peakE))/wfE);
                }
              }
              
  real[] deltG(real time,
              real[] peakG,
              real[] wrG,
              real[] wfG) {
                for (i in 1:size(t)){
                  return t<peakG ? exp((-(t-peakG)^2)/wrG) : exp((-(t-peakG))/wfG);
                }
              }
              
  real[] deltM(real time,
              real[] peakM,
              real[] wrM,
              real[] wfM) {
                for (i in 1:size(t)){
                  return t<peakM ? exp((-(t-peakM)^2)/wrM) : exp((-(t-peakM))/wfM);
                }
              }
  
  real[] TIV(real time,
            real[] y,  // system components {target cell, infected cell, viral load}
            real[] params, // parameters
            real[] x_r,
            int[] x_i) {
      
      real T = y[1];
      real I = y[2];
      real V = y[3];
      real N = x_i[1];
      
      real d = params[1];
      real delta = params[2];
      real c = params[3];
      real kt.Aa = params[4];
      real p = params[5];
      real ktauE = params[6];
      real ktauG = params[7];
      real ktauM = params[8];
      
      for (i in 1:t) {
      dT_dt = (params[1] * 1e+12) - (params[1] * y[1]) - (params[4] * y[1] * y[3]);
      dI_dt = (params[4] * y[1] * y[3]) - ((params[2] + params[6] * deltE[i]) * y[2]);
      dV_dt = (params[5] * y[2]) - ((params[3] + params[7] * deltG[i] + params[8] * deltM[i]) * y[3]); 
      
      return {dT_dt, dI_dt, dV_dt};
      }
    }

}

data {
  int<lower = 1> n_obs; // Number of obs
  real y0[3];
  real t0;
  int<lower = 1> n_params; // Number of model parameters
  int<lower = 1> n_difeq; // Number of differential equations in the system
  int<lower = 1> n_fake; // This is to generate "predicted"/"unsampled" data
  real ts[n_obs]; // Time points that were sampled
  int N; 
  real fake_ts[n_fake]; // Time points for "predicted"/"unsampled" data
}

transformed data {
  real x_r[0];
  int x_i[1] = {N};
}

parameters {
  real<lower = 0> params[n_params]; // Model parameters
  real<lower = 0> V0; // Initial viral load
}

transformed parameters{
  real y_hat[n_obs, n_difeq]; // Output from the ODE solver
  real y0[n_difeq]; // Initial conditions for both S and I
  {
  real params[8];
  real params[1] = d;
  real params[2] = delta;
  real params[3] = c;
  real params[4] = kt.Aa;
  real params[5] = p;
  real params[6] = ktauE;
  real params[7] = ktauG;
  real params[8] = ktauM;
  
  y_hat = integrate_ode_rk45(TIV, y0, t0, ts, params, x_r, x_i);
  }
}

model {
  //priors
  d ~ normal(0, 0.5); //truncated at 0
  delta ~ normal(0, 1);
  c ~ normal(0, 1);
  kt.Aa ~ normal(0, 1);
  p ~ normal(0, 1);
  ktauE ~ normal(0, 1);
  ktauG ~ normal(0, 1);
  ktauM ~ normal(0, 1);
  
  //sampling distribution
  //col(matrix x, int n) - The n-th column of matrix x. Here the number of infected people
  cases ~ neg_binomial_2(col(to_matrix(y), 2), phi);
}

generated quantities {
  // Generate predicted data over the whole time series:
  real fake_I[n_fake, n_difeq];
  
  fake_I = integrate_ode_rk45(SI, y0, t0, fake_ts, params, x_r, x_i);
  
}

Again, many thanks for any help offered!

EDIT: @maxbiostat edited this post for syntax highlighting.

1 Like
6

If you want to verify that the ODE is specified correctly in Stan, probably the easiest way is to use the expose_stan_functions function in rstan. That’ll allow you to call functions defined in Stan from R, and you can check them against your R implementation.

You’re limited to only using integrate_ode_rk45 though – I don’t think integrate_ode_bdf works with expose_stan_functions.

2 Likes
7

I’ll give that a shot, thank you!

8

As expected, I seem to be getting errors with defining the time-dependent parameters themselves. I need some kind of if/else structure, and I thought that using the format condition ? function if condition true : function if condition false would suffice, but instead R throws the error:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Binary infix operator <= with functional interpretation logical_lt requires arguments of primitive type (int or real), found left type=real[ ], right arg type=real.
No matches for:

  logical_lt(real[ ], real)

Available argument signatures for logical_lt:

  logical_lt(int, int)
  logical_lt(int, real)
  logical_lt(real, int)
  logical_lt(real, real)

However, according to Stan, time needs to be of type real[ ] .
My question is, does another if/else function exist in Stan that I can use?

9

Can you rewrite this comparison as a for loop?

for(n in 1:N) {
  if(array[n] < value) {
    ...
  } else {
    ...
  }
}
10

Yes, I apologize for the lack of clarity:

for(n in 1:length(time)) {
   if(time[n] < peak_time) {
   f <- exp((-((time[n]-peak_time)^2))/width_rising)
   } else {
   f <- exp((-(time[n]-peak_time))/width_falling)
  }
}

This comparison acts as a forcing function, where the authors can dictate the peak, rising, and falling width for the curves that represent the various components of the adaptive immune system.

I hope this makes my intentions clearer. Thank you for your continued help!

11

Did this answer your question here:

I think I missed that the first time through.

1 Like
12

Yes! So the code I responded with should do the trick?
Regardless, I have a follow-up question: are there any Stan packages that perform Approximate Bayesian Computation (ABC) using rejection sampling/SMC?

13

I’m not sure if it will solve the problem – that’s something you’ll have to check yourself. It just sounded to me like you were comparing an array to a scalar (something that works in R but not stan), and so I was just typing out the Stan way of doing that.

Another thing, I’m not totally sure about this, but since you have a forcing function with a discontinuous first derivative, you might need to split the ODE solve into two pieces. Like solve from t = 0 to t = peak_time, stop the solver, and then solve from t = peak_time onwards.

If you’re using the RK45 solver you can look up if non-smooth forcing functions are a problem (and if you’re using BDF or Adams look up the question for BDF solvers or Adams-Moulton solvers).

1 Like