In gmonette/Tcells: Fits Tcell dynamic models with mixed models

Life and death cell processes TODO TODO

Is this a link to Plans: link to Plans

Suppose that the population of a particular type of cell develops through a birth and death Poisson process in which each has probability $\beta$ or reproducing and a probability $\delta$ of dying in one unit of time. Assume that the unit is sufficiently small, hence also $\beta$ and $\delta$, that the two processes can be treated as independent. Let $X_0$ be the number of cells at time 0 and consider small increments of time $t_1$, $t_2$ and $t_3$. Let $d_1$, $d_2$, $d_3$ be the changes in the number of cells from time 0 to time $t_1$, from time $t_1$ to time $t_1 + t_2$ and from time $t_1 + t_2$ to time $t_1 + t_2 + t_3$ respectively.

Continue with "eqn limiting distribution of cell in ....." and formulas in photo (check consistency)

Model for T-cell dynamics

Subjects are observed at a number of occasions at different times. On each occasion the densities of cells of five different types are recorded.

The process for these five types, $N, S, C, T$ and $E$, are modelled as follows:

[\begin{aligned} {N_{t + 1}} = & {N_t} + {B_t} - N_t^d - N_t^t \ {S_{t + 1}} = & {S_t} + N_t^t - S_t^d - S_t^t \ {C_{t + 1}} = & {C_t} + S_t^t - C_t^d - C_t^t \ {T_{t + 1}} = & {T_t} + C_t^t - T_t^d - T_t^t \ {E_{t + 1}} = & {E_t} + T_t^t - E_t^d \ \end{aligned} ]

where

1. $X_t^d$ represents the number of cells of type $X$ that die between time $t$ and time $t+1$,
2. $X_t^t$ represents the number of cells of type $X$ that transition to the next type,
3. $B_t$ represents the number of new cells of type $N$.

Letting $d_t = \tau_{t+1} - \tau_t$ be the number of units of time elapsed between the observations at times $\tau_t$ and $\tau_{t+1}$, the hypothetical distribution of each component is:

[\begin{aligned} {B_t}\sim & Poisson(\lambda d_t) \ N_t^d\sim & Poisson\left( {{N_t}{\delta _N d_t}} \right)\ N_t^t\sim & Poisson\left( {{N_t}{\psi _N} d_t} \right)\ S_t^d\sim & Poisson\left( {{S_t}{\delta _S} d_t} \right) \ S_t^t\sim & Poisson\left( {{S_t}{\psi _S} d_t} \right) \ C_t^d\sim & Poisson\left( {{C_t}{\delta _C} d_t} \right) \ C_t^t\sim & Poisson\left( {{C_t}{\psi _C} d_t} \right) \ T_t^d\sim & Poisson\left( {{T_t}{\delta _T} d_t} \right)\ T_t^t\sim & Poisson\left( {{C_t}{\psi _T} d_t} \right) \ E_t^d\sim & Poisson\left( {{E_t}{\delta _E} d_t} \right) \ \end{aligned} ]

This model has ten parameters: $\lambda$ representing the expected number of new cells of type $N$ per unit of time, and nine parameters of the form $\delta_X, \psi_X$, representing the probability of death or transition per cell per unit of time.

If $d_t$ is almost constant, it can be be omitted in the model. For simplicity, it is omitted in the following formulas.

The values of some or all of these parameters are expected to vary from subject to subject. Using a mixed model in which the parameters are random allows this feature to be include in modeling.

The model should be accurate if the number of dying and transitioning cells of each type is relatively small compared with the total number of cells of that type within each observed time period. In this case, it may also be reasonable to assume that the $B_t, N_t^d, N_t^t, ...$ components are independent within occasion $t$ and sequentially conditionally independent, i.e. the distribution of $B_t |N_t$ is independent of the distribution of $B_{t+1}|N_{t+1}$, etc.

The change in cell counts between occasions is: [\Delta \left[ {\begin{array}{{20}{c}} {{N_{t + 1}}} \ {{S_{t + 1}}} \ {{C_{t + 1}}} \ {{T_{t + 1}}} \ {{E_{t + 1}}} \end{array}} \right] = A \left[ {\begin{array}{{20}{c}} {{B_t}} \ {N_t^d} \ {N_t^t} \ {S_t^d} \ {S_t^t} \ {C_t^d} \ {C_t^t} \ {T_t^d} \ {T_t^t} \ {E_t^d} \end{array}} \right]] where [A = \left[ {\begin{array}{{20}{c}} 1&{ - 1}&{ - 1}&0&0&0&0&0&0&0 \ 0&0&1&{ - 1}&{ - 1}&0&0&0&0&0 \ 0&0&0&0&1&{ - 1}&{ - 1}&0&0&0 \ 0&0&0&0&0&0&1&{ - 1}&{ - 1}&0 \ 0&0&0&0&0&0&0&0&1&{ - 1} \end{array}} \right]] The expectation of the change in cell counts is: [E\left( {\Delta \left[ {\begin{array}{{20}{c}} {{N_{t + 1}}} \ {{S_{t + 1}}} \ {{C_{t + 1}}} \ {{T_{t + 1}}} \ {{E_{t + 1}}} \end{array}} \right]} \right) = \left[ {\begin{array}{{20}{c}} {\lambda - {\delta _N}{N_t} - {\psi _N}{N_t}} \ {{\psi _N}{N_t} - {\delta _S}S - {\psi _S}{S_t}} \ {{\psi _S}{S_t} - {\delta _C}{C_t} - {\psi _C}{C_t}} \ {{\psi _C}{C_t} - {\delta _T}{T_t} - {\psi _T}{T_t}} \ {{\psi _T}{T_t} - {\delta _E}{E_t}} \end{array}} \right] = \left[ {\begin{array}{{20}{c}} {\lambda - {N_t}\left( {{\delta _N} + {\psi _N}} \right)} \ {{\psi _N}{N_t} - {S_t}\left( {{\delta _S} + {\psi _S}} \right)} \ {{\psi _S}{S_t} - {C_t}\left( {{\delta _C} + {\psi _C}} \right)} \ {{\psi _C}{C_t} - {T_t}\left( {{\delta _T} + {\psi _T}} \right)} \ {{\psi _T}{T_t} - {\delta _E}{E_t}} \end{array}} \right] ] and the variance-covariance matrix is [\operatorname{Var} \left( {\Delta \left[ {\begin{array}{{20}{c}} {{N_{t + 1}}} \ {{S_{t + 1}}} \ {{C_{t + 1}}} \ {{T_{t + 1}}} \ {{E_{t + 1}}} \end{array}} \right]} \right) = A\left[ {diag\left( {\begin{array}{{20}{c}} \lambda &{{\delta _N}{N_t}}&{{\psi _N}{N_t}}&{{\delta _S}{S_t}}&{{\psi _S}{S_t}}&{{\delta _C}{C_t}}&{{\psi _C}{C_t}}&{{\delta _T}{T_t}}&{{\psi _T}{T_t}}&{{\delta _E}{E_t}} \end{array}} \right)} \right]A']

[ = \left[ {\begin{array}{*{20}{c}} {\lambda + {N_t}\left( {{\delta _N} + {\psi _N}} \right)}&{ - {\psi _N}{N_t}}&0&0&0 \ { - {\psi _N}{N_t}}&{{\psi _N}{N_t} + {S_t}\left( {{\delta _S} + {\psi _S}} \right)}&{ - {\psi _S}{S_t}}&0&0 \ 0&{ - {\psi _S}{S_t}}&{{\psi _S}{S_t} + {C_t}\left( {{\delta _C} + {\psi _C}} \right)}&{ - {\psi _C}{C_t}}&0 \ 0&0&{ - {\psi _C}{C_t}}&{{\psi _C}{C_t} + {T_t}\left( {{\delta _T} + {\psi _T}} \right)}&{ - {\psi _T}{T_t}} \ 0&0&0&{ - {\psi _T}{T_t}}&{{\psi _T}{T_t} + {\delta _E}{E_t}} \end{array}} \right]]

This formulation of the model allows an unconstrained parametrization by using the log of each parameter, each of which is positive, as the actual working parameter for optimization.

Equilibrium and stability conditions

[E\left( {\Delta \left[ {\begin{array}{{20}{c}} {{N_{t + 1}}} \ {{S_{t + 1}}} \ {{C_{t + 1}}} \ {{T_{t + 1}}} \ {{E_{t + 1}}} \end{array}} \right]} \right) = \left[ {\begin{array}{{20}{c}} {\lambda - {\delta _N}{N_t} - {\psi _N}{N_t}} \ {{\psi _N}{N_t} - {\delta _S}S - {\psi _S}{S_t}} \ {{\psi _S}{S_t} - {\delta _C}{C_t} - {\psi _C}{C_t}} \ {{\psi _C}{C_t} - {\delta _T}{T_t} - {\psi _T}{T_t}} \ {{\psi _T}{T_t} - {\delta _E}{E_t}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\lambda - {N_t}\left( {{\delta _N} + {\psi _N}} \right)} \ {{\psi _N}{N_t} - {S_t}\left( {{\delta _S} + {\psi _S}} \right)} \ {{\psi _S}{S_t} - {C_t}\left( {{\delta _C} + {\psi _C}} \right)} \ {{\psi _C}{C_t} - {T_t}\left( {{\delta _T} + {\psi _T}} \right)} \ {{\psi _T}{T_t} - {\delta _E}{E_t}} \end{array}} \right] ]

[E\left( {\Delta \left[ {\begin{array}{{20}{c}} {{N_{t + 1}}} \ {{S_{t + 1}}} \ {{C_{t + 1}}} \ {{T_{t + 1}}} \ {{E_{t + 1}}} \end{array}} \right]} \right) = \left[ {\begin{array}{{20}{c}} {\lambda - {N_t}\left( {{\delta _N} + {\psi _N}} \right)} \ {{\psi _N}{N_t} - {S_t}\left( {{\delta _S} + {\psi _S}} \right)} \ {{\psi _S}{S_t} - {C_t}\left( {{\delta _C} + {\psi _C}} \right)} \ {{\psi _C}{C_t} - {T_t}\left( {{\delta _T} + {\psi _T}} \right)} \ {{\psi _T}{T_t} - {\delta _E}{E_t}} \end{array}} \right] = \left[ {\begin{array}{{20}{c}} \lambda &{ - \left( {{\delta _N} + {\psi _N}} \right)}&0&0&0&0 \ 0&{{\psi _N}}&{ - \left( {{\delta _S} + {\psi _S}} \right)}&0&0&0 \ 0&0&{{\psi _S}}&{ - \left( {{\delta _C} + {\psi _C}} \right)}&0&0 \ 0&0&0&{{\psi _C}}&{ - \left( {{\delta _T} + {\psi _T}} \right)}&0 \ 0&0&0&0&{{\psi _T}}&{ - {\delta _E}} \end{array}} \right]\left[ {\begin{array}{{20}{c}} 1 \ {{N_t}} \ {{S_t}} \ {{C_t}} \ {{T_t}} \ {{E_t}} \end{array}} \right]]

Adjustment for cell concentrations from blood samples of different volumes

If the response variable is either a number of cells counted in a situation in which the volume of the sample from which cells are counted varies from observation to observation or if the response is a cell concentration, i.e. number of cells observed divided by the volume of the sample examined, then the model above needs to be adjusted.

Let $X_t, Y_t$ be cells concentrations observed in a smaple of volume $V_t$. Let $X_t^c, Y_t^c$ be the actual counts so that [\left( {\begin{array}{{20}{c}} {{X_t}} \ {{Y_t}} \end{array}} \right) = \frac{1}{{{V_t}}}\left( {\begin{array}{{20}{c}} {X_t^c} \ {Y_t^c} \end{array}} \right)] If $\left( {\begin{array}{{20}{c}} {X_t^c} \ {Y_t^c} \end{array}} \right)$ has expectation $\left( {\begin{array}{{20}{c}} {{\mu X}} \ {{\mu _Y}} \end{array}} \right)$ and variance $\left[ {\begin{array}{{20}{c}} {{\sigma {XX}}}&{{\sigma {XY}}} \ {{\sigma {YX}}}&{{\sigma {YY}}} \end{array}} \right]$ then the expectation and variance of $(X_t, Y_t)'$ will be $\frac{1}{{{V_t}}}\left( {\begin{array}{{20}{c}} {{\mu _X}} \ {{\mu _Y}} \end{array}} \right)$ and $\frac{1}{{V_t^2}}\left[ {\begin{array}{*{20}{c}} {{\sigma {XX}}}&{{\sigma {XY}}} \ {{\sigma {YX}}}&{{\sigma _{YY}}} \end{array}} \right]$, respectively.

Thus, blood sample volumes that vary between measurements can be taken into account but need to be incorporated into the model. Simply standardizing by, for example, multiplying the concentrations so they refer to a common volume will not work since the resulting relationship between the expectation and the variance of the response variables will not be that obtained above assuming Poisson components.

If the measurements are cell concentrations that are always obtained from the same blood volume, then the easiest approch is to multiply each obervation by the common value of $V$ which will yield the cell counts modeled as discussed in the previous sections.

Using mixed models

In most applications of mixed models, there is a univariate response variable whose expectation is determined by a linear model some of whose coefficients are fixed unknown parameters and others vary randomly from cluster to cluster. The model for the variance of the random parameters can be selected from a limited number of possibilies. Within each cluster, responses are independent with the same conditional within cluster variance for all observations.

This application is particularly interesting in that none of these common assumptions are reasonable.

The response is multivariate (the five cell types) and the variance is different at each occasion depending on parameters of the expectation model.

Few packages have the flexibility to accomodate these issues. The 'nlme' package has facilities to deal with:

1. Different models for the variance of random parameters: However it does not include models that would be reasonable for this kind of problem. A full free variance model would have too many parameters so that some independence assumptions are necessary. However, the only such model readily available in 'nlme' forces independence from the intercept random effect which imposes a very unrealistic assumption on the model.
2. Variances that vary depending on a fixed covariate or on the expected value of the response: However, there is no provision for variance that depends in a complex way on parameters of the linear model that are estimated iteratively in fitting the model.
3. Correlations between responses: However, again, there is no provision for correlations that depend on linear parameters nor for the type of correlation pattern in this problem.

The virtue of 'nlme' is that its source is open and it is possible to create new 'methods' to implement models for the random parameter variance, the conditional response variance and the within-occasion response covariances that occur in this model.

The first part, which was more 'theoretically' challenging is complete. The other two parts are well under way.

The software to implement this is being built in a R package called 'Tcells' on github. It can be installed in R with the following commands:

install.packages('devtools')
library(devtools)
install_github('gmonette/Tcells')
library(Tcells)

Since the package is being developed it is good practice to reinstall it frequently.

gmonette/Tcells documentation built on May 17, 2019, 7:25 a.m.