Age structured model to track cell-age and ki67 exprerssion in naive T cells
11 March, 2020
PDE Model
In order to track the dynamics of cells of age ‘a’ expressing the nuclear protein Ki67 with an intensity ‘k’ through time `t’, we use the following PDE, \[\begin{equation}
\frac{\partial u}{\partial t} +\frac{\partial u}{\partial a} - \beta \, k \, \frac{\partial u}{\partial k} = - \delta(a) \, u.
\label{eq:PDE}
\end{equation}\] - \(\beta\) is the rate of loss of Ki67 expression, and
- \(\delta\) is the rate with which cells are lost from the population either by death or by differentiation. We assume that it varies with cell-age.
To minimise the experimental variation we normalise the intensity of Ki67 such that \(k\) \(\in (0, 1)\). The Ki67 intensity is highest (\(k=1\)) when cells enter division and then decays exponentially with rate constant \(\beta\). The threshold Ki67 intensity \(\bar{k}\) identifies the cells that are gated as Ki67\(^{+}\) in experimental observations (k \(\ge\) \(\bar{k}\)).
To solve this linear first-order PDE (eq.), we need to transform it into an ODE such as, \[\begin{equation} \frac{d}{ds}u(t(s), a(s), k(s)) = F(u, t(s), a(s), k(s)), \label{eq:ODE} \end{equation}\] along the characteristic curve (t(s), a(s), k(s)).
The solution of the ODE (eq.~) is given by, \[\begin{eqnarray} \begin{aligned} &\frac{du}{ds} = - \delta(a) \, u, \\ &\text{which can be solved to give,} \\ &u(s) = u_0 \, exp(- \int \delta(s) ds). \end{aligned} \label{eq:SOL} \end{eqnarray}\] Therefore, if we know \(u(t_0, a_0, k_0) \equiv \, u_0\) we can find \(u(t(s), a(s), k(s)\) using eq. , since (t(s), a(s), k(s)) and (t(0), a(0), k(0)) both lie on the same characteristic curve.
We track the folllowing three cohorts of naive T cells that form the boundary conditions for the PDE in eq. alomg this charcteristic curve.
- Initial population present at \(t_0 \Rightarrow u_\text{init} = u(t=0, a, k)\),
- Cells coming from the source \(\Rightarrow u_\text{source} = u(t, a=0, k)\) and
- Cells that entered division \(\Rightarrow u_\text{div} = u(t, a, k=1)\).
Total size of the population at time t is derived by integrating the sum of cell densities of these three cohorts over all allowable cell ages and Ki67 intensities. \[\begin{equation}
u(t) = \int_{a = 0}^t \int_{k = 0}^1 u_\text{init} + u_\text{source} + u_\text{div} \, \, dk \, \, da
\end{equation}\]
The proportions of Ki67\(^{+}\) cells within the naive compartment at time t is estimated by integrating the sum of cell densities of the three cohorts from the threshold Ki67 intensity \(\bar{k}\) to 1 and over all allowable cell ages. \[\begin{equation}
\kappa(t) = \frac{ \int_{a = 0}^t \int_{k = \bar{k}}^1 u_\text{init} + u_\text{source} + u_\text{div} \, \, dk \, \, da.}{u(t)}
\end{equation}\]
Thymic export
Mature CD4 SP cells migrate to the respective peripheral naive compartments continuously. The numbers of new naive T cells exported to the peripheral lymphoid organs varies across the mouse life-span, due to the changes in the SP population sizes owing to thymic involution. Moreover, the proportion of Ki67\(^{+}\) cells within the thymic SP CD4 compartment changes with host-age. In order to characterise the influx of Ki67\(^{+}\) and Ki67\(^{-}\) thymic emigrants into the peripheral naive compartments we define phenomenological function to track the counts and fraction of Ki67\(^{+}\) cells within the thymic SP CD4 subset.
Total counts of the thymic CD4 SP subset
The total numbers of CD4 SP cells increase rapidly upto 6-7 weeks of age and then drop gradually over time. We use a descriptor function (\(\theta(t)\), eq.~) to capture the dynamics of thymic SP cells varying with mouse age. \[\begin{equation} \theta(t) = \theta_{0} + (\theta_{f} \, t^{n}) \, \bigg(1 - \frac{t^{q}}{X^{q} + t^{q}}\bigg). \label{eq:Source_counts} \end{equation}\] The parameters \(\theta_{0}, \theta_{f}, X, n\) and \(q\) are estimated from the fits to the total counts of thymic CD4 and CD8 SP cells, separatley.
Dynamics of Ki67 expression within the thymic CD4 SP subsest
We fit an exponential function \(\epsilon\) that describes the time-dependent changes in the proportions Ki67\(^{+}\) cells (eq.~) to the fractions of Ki67\(^{+}\) and Ki67\(^{-}\) cells within the thymic efflux. \[\begin{equation} \epsilon(t) = \epsilon_{0} + e^{\epsilon_{f} \, (t + A)}. \label{eq:Source_ki67} \end{equation}\] The parameters \(\epsilon_{0}, \epsilon_{f}\) and \(A\) are estimated from the fits to the proportions of Ki67\(^{+}\) thymic CD4 SP cells.
# the function for thymic SP population -- changes with time
theta_t <- function(dpt){
## parameters estimated from spline fit to the timecourse of counts of source compartment -- SP
if (Population == "CD4") {
spline_est_list <- list("theta_0" = 4.3e5, 'theta_f' = 1.8e3, "n" = 2.1, "X" = 30, "q" = 3.7)
} else if (Population == "CD8") {
spline_est_list <- list("theta_0" = 9e4, 'theta_f' = 68, "n" = 3, "X" = 25, "q" = 4.25)
}
## days post t_0 -- dpt
## change in SP counts with time
theta = spline_est_list$theta_0 + (spline_est_list$theta_f * dpt^spline_est_list$n) *
(1 - ((dpt^spline_est_list$q)/((spline_est_list$X^spline_est_list$q) + (dpt^spline_est_list$q))))
return(theta)
}
# function for ki67 proportions in thymic output -- changes with time
eps_func <- function(Time){
## spline fitted separately to the proportions of ki67+ cells within thymic SP4 cells
## parameters estimated from spline fit to the timecourse of ki67 proportions within the source compartment -- SP
if (Population == "CD4") {
spline_est_list <- list("eps_0" = 0.13732103, 'eps_f' = 0.03337899, "A" = 2.92554110)
} else if (Population == "CD8") {
spline_est_list <- list("eps_0" = 0.24510453, 'eps_f' = 0.01559996, "A" = 14.83715328)
}
eps = exp(-spline_est_list$eps_f * (Time + spline_est_list$A)) + spline_est_list$eps_0
return(eps)
}
p1 <- ggplot()+
geom_line(aes(ts_pred, theta_t(ts_pred)), size =1.5, col = "#D87534") +
labs(x = 'Host age', y = NULL, title = paste0("Total ", Population , " SP pool size")) +
scale_y_log10(labels = fancy_scientific) + myTheme
p2 <- ggplot()+
geom_line(aes(ts_pred, eps_func(ts_pred)), size =1.5, col = "#D87534") +
labs(x = 'Host age', y = NULL, title = paste0("The proportion of Ki67high ", Population, " SP cells")) + myTheme
gridExtra::grid.arrange(p1, p2, nrow = 1)
Daily influx from thymic CD4 SP subset
We assume that the proportion of CD4 SP T cells entering naive pool daily, defined as \(\psi\), remains constant throughout the mouse life. We also consider that both the Ki67\(^{+}\) and Ki67\(^{-}\) SP cells follow the same kinetics of migration.
Therefore, the number of cells entering the total, Ki67\(^{+}\) and Ki67\(^{-}\) naive T cell subsets are given by, \[\begin{eqnarray} \begin{aligned} &\phi(t) = \psi \, \theta(t), \\ &\phi^{+}(t) = \psi \, \theta(t) \, \epsilon(t) \, \text{ and} \\ &\phi^{-}(t) = \psi \, \theta(t) \, (1 - \epsilon(t)) \, \text{, respectively. } \end{aligned} \label{eq:Source_influx} \end{eqnarray}\]
The distribution of Ki67 intensities within the thymic export
We track the distribution of Ki67 intensities within the thymic output using the \(\epsilon\) function that described the dynamics of Ki67 proportions in the thymus. \[\begin{equation} \kappa_{\theta}(k, t) = \frac{\epsilon(t)}{1 - \bar{k}} + \frac{\frac{1 - \epsilon(t)}{\bar{k}} - \frac{\epsilon(t)}{1 - \bar{k}}}{1+ \big(\frac{k}{\bar{k}}\big)^n} \label{eq:ki67_dist} \end{equation}\] By setting a high value for the Hill-exponent \((n = 40)\), we approximate the distribution in eq.~ to integrate to \(\epsilon(t)\) across \(\bar{k}\) to 1. Scaling the \(\kappa_{\theta}\) by multiplying it with the total number of cells that are exported from the thymus, we can estimate the numbers of Ki67\(^{+}\) cells withing the source influx. \[\begin{equation} \phi_{\kappa}(k, t) = \phi(t) \, \kappa_{\theta}(k, t) \label{eq:phi_ki67} \end{equation}\] Integrating \(\phi_{\kappa}(k, t)\) from 0 to 1 gives total thymic output (\(\phi(t)\)) and from \(\bar{k}\) to 1 gives the proprtion of Ki67\(^{+}\) cells within the thymic output i.e. \(\phi^{+}(t)\).
In our system, we define \(\bar{k}\) as the thresold Ki67 intensity used to gate as Ki67\(^{+}\) cells. We assume that the intensity of ki67 drops with a constant rate \(\beta\) with time \(k = e^{\beta \, t}\) (for normalised k, where \(k_0 = 1\)). Therefore, at time \(= \frac{1}{\beta} \rightarrow k = \bar{k} \text{ and } \rightarrow \bar{k} = \frac{1}{e}\).
## the function for thymic influx into the naive T cell compartment
phi_t <- function(Time, parms) {
N0 = parms[1]
#q1 = parms[2]
t0 = 5
## psi is the proportion of SP cells that enter naive pool daily -- constant here
#psi = (exp(N0) * q1)/(theta_t(t_0) * (exp(q1 * t_0) - 1))
### simpler assumption is that intital age and ki67 distribution is same as the thymic output -- because t0 = 5 days.
### g(a, k) = phi(t0, k) --> integral of g(a, k) across 0 to t0 == 5 * phi_0 == 5 * theta_0 * psi.
### therefore psi = N0/(5 * theta_0) --- this wouldnt hold if t0 was even slightly larger.
psi = exp(N0)/(5 * theta_t(t0))
value = ifelse(Time >= t0, psi * theta_t(Time),
0)
## the proportional influx at each time 't' is,
return(value)
}
# distribution of ki67 expression in thymic output -- changes with time
k_dist <- function(ki, Time){
k_bar = 1/exp(1)
n = 60
t_0 = 5
value = (eps_func(Time)/(1 - k_bar)) +
((((1 - eps_func(Time))/(k_bar)) - (eps_func(Time)/(1 - k_bar)))/(1 + (ki/k_bar)^n))
return(value)
}
## function to give ki67 distribution of thymic output at a given time 't' -- change with time.
phi_k_t <- function(ki, Time, parms){
t0 = 5
value = phi_t(Time, parms) *
k_dist(ki, Time)
return(value)
}
## integrating phi_k_t from 0 to 1 gives total thymic output
# from k_bar (threshold for ki67+ signal) to 1 gives the total numbers of ki67+ cells in the thymic output
phi_int <- function(Time, parms, low_lim_kbar){
value <- function(Time, parms, low_lim_kbar) {
k_bar = 1/exp(1)
low_lim = ifelse(low_lim_kbar == FALSE, 0,
k_bar)
return(integrate(phi_k_t, lower = low_lim, upper = 1, Time=Time, parms=parms)$value)
}
return(mapply(value, Time, MoreArgs = list(parms, low_lim_kbar)))
}
phi_df <- data.frame("Time" = ts_pred,
"Total" = phi_int(ts_pred, start_vec, F),
"Ki_pos" = phi_int(ts_pred, start_vec, T)) %>%
mutate(Ki_neg = Total - Ki_pos) %>%
gather(-Time, key = Subpopulation, value = Counts)
ggplot(phi_df) +
geom_line(aes(Time, Counts/1e6, col = Subpopulation), size =1.5) +
#scale_y_log10(labels = fancy_scientific) +
scale_color_discrete(breaks = c("Total", "Ki_pos", "Ki_neg"),
labels = c( "Total", "Ki_pos" = "Ki67high", "Ki_neg" = "Ki67low")) +
labs(x = 'Host age', y = NULL, title = "Thymic output (numbers in millions)") + myTheme
Tracking the cohorts of cells of age ‘a’ and Ki67 intensity ‘k’
The initial population
We define the distribution of cell-ages and ki67 intensities of naive T cells pre-existing at time t0 using a function \(g(a, k)\). \[\begin{equation}
u_\text{init}(t, a, k) = g(a, k).
\end{equation}\]
The total size of this initial population is given by the integral of their age and ki67 distribution across all possible values of cell ages (0 to \(t_0\)) and ki67 intensities (0 to 1). Similarly, the proportion of Ki67\(^{+}\) cells at time \(t_0\) is given by dividing the integral of \(g(a, k)\) across all ages and from \(\bar{k}\) to 1 with the total pool size \(N_0\). \[\begin{eqnarray} N_0 = \int_0^{t_0} \int_0^1 g(a, k) \, dk \, da \\ \kappa_0 = \frac{\int_0^{t_0} \int_{\bar{k}}^1 g(a, k) \, dk \, da}{N_0} \end{eqnarray}\]
The initial population of naive T cells includes the cohort of thymic emigrants that entered the periphery at \(t_0, \, \rightarrow g(0, k) = \phi(t_0, k).\) Since, the earliest time-point of experimental observations is day 5 since birth, we assume that the effects of cell-age on naive T cells dynamics would be extremely small and therefore, consider that the \(g(a, k)\) is uniform and more or less reflects the distributions in the thymic output i.e., \(\phi(t_0, k)\).
Expressing the rate of daily thymic influx in terms of pool size of initial population – \(N_0\).
\[\begin{eqnarray} \begin{aligned} &N_0 = \int_{s=0}^{t_0} g(s) \, ds \quad = \int_{s=0}^{t_0} \phi(t_0) \, ds \quad = \psi \, \theta(t_0) \, t_0 \\ &\therefore \psi = \frac{N_0}{t_0 \, \theta(t_0)}. \end{aligned} \end{eqnarray}\]
## tracking preexisting cells that were present at t0
## the age distribution of cells presnt at t0
g_age <- function(age, parms){
q1 = 0 # parms[2] ## We asssumed it to be flat since t0 is very small.
#### simpler assumption that age and ki67 distr. are flat because t0 is very small
t_0 = 5
# g(a)
value = ifelse(age >= 0,
phi_t(t_0, parms) * exp(q1 * age),
0)
return(value)
}
## the age and ki67 distribution of cells presnt at t0
g_k_a <- function(ki, age, Time, parms){
# g(a, k)
value = g_age(age, parms) * k_dist(ki, Time) ### proportional reduction in cells of age 'a' and ki67 'k'
return(value)
}Along the charcteristic curve \((a(s), t(s), k(s))\) the PDE is reduced to an ODE. By solving the charactristics system of ODEs the we determine the solution along the first boundary condition i.e. the initial population as follows. \[\begin{eqnarray} \begin{aligned} &\frac{dt}{ds} = 1 \quad \rightarrow t = s + t_0 \quad \rightarrow t_0 = 0 \quad \rightarrow s = t \\ &\frac{da}{ds} = 1 \quad \rightarrow a = s + a_0 \quad \rightarrow a_0 = a - t \\ &\frac{dk}{ds} = -\beta \, k \quad \rightarrow k = k_0 \, e^{-\beta \, s} \quad \rightarrow k_0 = k \, e^{\beta \, t}. \end{aligned} \end{eqnarray}\] Therefore, \[\begin{equation} u_\text{init}(0, a_0, k_0) = \phi(0, k \, e^{\beta \, t}). \end{equation}\]
Cellular Turnover: We assumed that the rate wth which cells are lost either by death or differentiation i.e., ‘\(\delta\)’ changes exponentially with their age, \[\begin{equation} \delta (a) = \delta_0 \, e^{-r \, a}. \end{equation}\]
### cellular turnover changing with cellage
delta_age <- function(age, parms){
delta0 = parms[2]
r_del = parms[3]
return(delta0 * exp(-r_del * age))
}
delta_integ <- function(lo_lim, up_lim, parms){
delta_value <- function(lo_lim, up_lim, parms){
integrate(delta_age, lower = lo_lim, upper = up_lim, parms=parms)$value
}
## mapply return the out put of the value function across each elements of vectors
## lower and upper for the same values of params as a list
return( mapply(delta_value, lo_lim, up_lim, MoreArgs = list(parms)))
}Using \(u_{p}(t_0, a_0, k_0)\) and \(\delta(a)\) the initial population can be tracked on the characteristic curve `\(S\)’, \[\begin{equation} \begin{aligned} u_\text{init}(t, a, k) &= u_\text{init}(0, a_0, k_0) \, exp\big(- \int_{a-t}^{a} \delta(\tau) \; d\tau \big)\\ &= \phi\big(0, k \, e^{\beta \, t}\big) \, exp\big(- \int_{a-t}^{a} \delta(\tau) \; d\tau \big) \end{aligned} \end{equation}\]
The distribution of ki67 expression of cells within the initial popualtion at time t is, \[\begin{equation}
u_\text{init}(t, k) = \int_0^t \phi\big(0, k \, e^{\beta \, t}\big) \, exp\big(- \int_0^{t} \delta(\tau + a - t) \; d\tau\big) \; da.
\end{equation}\]
Therefore the total size (\(N_\text{p}\)) of the ‘Pre-existing cohort’ and the porportion of Ki67\(^+\) cells (\(\kappa_\text{p}\)) within it are given by, \[\begin{eqnarray} \begin{aligned} & N_\text{init}(t) = \int_0^1 u_\text{p}(t, k) \; dk \\ & \kappa_\text{init}(t) = \frac{\int_{\bar{k}}^1 u_\text{p}(t, k) \, dk }{N_\text{p}(t)}. \end{aligned} \end{eqnarray}\]
### pre-existing cells age and ki67 distribution
N_init_k_a <- function(age, ki, Time, parms){
delta0_log = parms[2] ### log transforming parameter ---> delta0 > 0
delta0 = exp(delta0_log)
r_del = parms[3] ### free parameter
t0 = 5
#value = g_k_a(ki, age,Time, parms) * exp(- delta_integ(lo_lim = 0, up_lim = age, parms))
### analytically solved inetgral for delta_age to speed up the estimation
value = g_k_a(ki, age - Time + t0, Time, parms) * exp(- (delta0/r_del) * (exp(-r_del * (age - Time + t0))
- exp(-r_del * age)))
return(value)
}
## ki67 distribution of the cohorts entered from the thymus post 't0'
## integral across all possible ages
N_init_k <- function(ki, Time, parms){
init_value <- function(ki, Time, parms){
t0 = 5
return(integrate(N_init_k_a, lower = Time - t0, upper = Time, ki=ki, Time=Time, parms = parms)$value)
}
## mapply return the out put of the value function across each elements of vectors
## k and t for the same values of params as a list
return(mapply(init_value, ki, Time, MoreArgs = list(parms)))
}
## number of ki67+ cells within the cohorts entered from the thymus post 't0' at time 't'
k_init <- function(Time, parms){
k_value <- function(Time, parms){
k_bar = 1/exp(1)
## inetgrating across k_bar to 1
return(integrate(N_init_k, lower = k_bar, upper = 1, Time=Time, parms=parms)$value)
}
## mapply return the out put of the value function across each elements of vector
## t for the same values of params as a list
return(mapply(k_value, Time, MoreArgs = list(parms)))
}
## total population size of the cohorts entered from the thymus post 't0' at time 't'
N_init <- function(Time, parms){
init_value <- function(Time, parms){
## inetgrating across all k values
return(integrate(N_init_k, lower = 0, upper = 1, Time=Time, parms=parms)$value)
}
## mapply return the out put of the value function across each elements of vector
##t for the same values of params as a list
return(mapply(init_value, Time, MoreArgs = list(parms)))
}
## proportions of ki67+ cells within the cohorts entered from the thymus post 't0' at time 't'
ki_prop_init <- function(Time, parms){
# numbe rof ki67+ cells/ totla population size
value =
k_init(Time, parms)/
N_init(Time, parms)
return(value)
}# integration
time_integ_init <- paste0("time taken for integration -- ",
system.time(
c(counts_init_vec <- mcmapply(N_init, ts_pred, MoreArgs =list(start_vec), mc.cores = num_cores),
ki_prop_init_vec <- mcmapply(ki_prop_init, ts_pred, MoreArgs =list(start_vec), mc.cores = num_cores)),
gcFirst = TRUE)[3]
)
init_res_df <- data.frame("time" = ts_pred,
"Counts" = counts_init_vec,
"Ki67high proportions" = ki_prop_init_vec) %>%
gather(-time, key = "features", value = "results")
ggplot()+
geom_line(data = init_res_df, aes(time, results), col = "#D87534", size = 1.5) +
labs(title = paste0("Initial cohort - naive ", Population, " T cells"), x = "Host age (days)", y = NULL) +
facet_wrap( .~ features, scales = "free") + myTheme
print(time_integ_init)Cells coming from the source
We assume that cells of age 0 constantly enter the target population from the source compartment which sets the second boundary condition for the PDE in eq. , \[\begin{equation*}
u_\text{source}(t, a=0, k) = \phi(t, k).
\end{equation*}\] The characteristics system of ODEs for this boundary cod=ndition is, \[\begin{eqnarray}
\begin{aligned}
&\frac{da}{ds} = 1 \quad \rightarrow a = s + a_0 \quad \rightarrow a_0 = 0 \quad \rightarrow s = a \\
&\frac{dt}{ds} = 1 \quad \rightarrow t = s + t_0 \quad \rightarrow t_0 = t - a \\
&\frac{dk}{ds} = -\beta \, k \quad \rightarrow k = k_0 \, e^{-\beta \, s} \quad \rightarrow k_0 = k \, e^{\beta \, a}.
\end{aligned}
\end{eqnarray}\] Therefore, \[\begin{equation}
u_\text{source}(t_0, 0, k_0) = \phi(t_0, k_0)= \phi\big(t-a, k \, e^{\beta \, a}\big).
\end{equation}\]
Using \(\phi(t_0,k_0)\) and \(\delta(a)\) the cohort of cells (age = 0) entering from the source can be tracked across time t, \[\begin{equation}
\begin{aligned}
u_\text{source}(t, a, k) &= u_\text{source}(t_0, 0, k_0) \, exp\big(- \int_0^a \delta(\tau) \; d\tau \big)\\
&= \phi\big(t-a, k \, e^{\beta \, a}\big) \, exp\big(- \int_0^a \delta(\tau) \; d\tau \big)
\end{aligned}
\end{equation}\] The distribution of ki67 expression of the cohorts entered from the thymus post ‘\(t_0\)’ at time t is, \[\begin{equation}
u_\text{source}(t, k) = \int_0^t \phi\big(t-a, k \, e^{\beta \, a}\big) \, exp\big(- \int_0^a \delta(\tau) \; d\tau\big) \; da.
\end{equation}\]
Therefore the total size (\(N_\text{source}\)) of the cohort of cells derived from the source compartment and the porportion of Ki67\(^+\) cells (\(\kappa_\text{source}\)) within them are given by, \[\begin{eqnarray} \begin{aligned} & N_\text{source}(t) = \int_0^1 u_\text{s}(t, k) \; dk \\ & \kappa_\text{source}(t) = \frac{\int_{\bar{k}}^1 u_\text{s}(t, k) \, dk }{N_\text{s}(t)}. \end{aligned} \end{eqnarray}\]
### cells coming from source -- age 0
N_theta_k_a <- function(age, ki, Time, parms){
delta0_log = parms[2] ### log transforming parameter ---> delta0 > 0
delta0 = exp(delta0_log)
r_del = parms[3] ### free parameter
t0 = 5
#value = phi_k_t(ki * exp(beta * age), Time - age, parms) * exp(- delta_integ(lo_lim = 0, up_lim = age, parms))
### analytically solved inetgral for delta_age to speed up the estimation
value = phi_k_t(ki * exp(beta * age), Time - age, parms) * exp(- (delta0/r_del) * (1 - exp(- r_del * age)))
return(value)
}
## ki67 distribution of the cohorts entered from the thymus post 't0'
N_theta_k <- function(ki, Time, parms){
theta_value <- function(ki, Time, parms){
t0 =5
return(integrate(N_theta_k_a, lower = 0, upper = Time - t0, ki=ki, Time=Time, parms = parms)$value)
}
## mapply return the out put of the value function across each elements of vectors
## k and t for the same values of params as a list
return(mapply(theta_value, ki, Time, MoreArgs = list(parms)))
}
## number of ki67+ cells within the cohorts entered from the thymus post 't0' at time 't'
k_theta <- function(Time, parms){
k_value <- function(Time, parms){
k_bar = 1/exp(1)
## inetgrating across k_bar to 1
return(integrate(N_theta_k, lower = k_bar, upper = 1, Time=Time, parms=parms)$value)
}
## mapply return the out put of the value function across each elements of vector
## t for the same values of params as a list
return(mapply(k_value, Time, MoreArgs = list(parms)))
}
## total population size of the cohorts entered from the thymus post 't0' at time 't'
N_theta <- function(Time, parms){
theta_value <- function(Time, parms){
## inetgrating across all k values
return(integrate(N_theta_k, lower = 0, upper = 1, Time=Time, parms=parms)$value)
}
## mapply return the out put of the value function across each elements of vector
##t for the same values of params as a list
return(mapply(theta_value, Time, MoreArgs = list(parms)))
}
## proportions of ki67+ cells within the cohorts entered from the thymus post 't0' at time 't'
ki_prop_theta <- function(Time, parms){
# numbe rof ki67+ cells/ totla population size
value =
k_theta(Time, parms)/
N_theta(Time, parms)
return(value)
}#Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
# integration
time_integ_source <- paste0("time taken for integration -- ",
system.time(
c(counts_theta_vec <- mcmapply(N_theta, ts_pred, MoreArgs =list(start_vec), mc.cores = num_cores),
ki_prop_theta_vec <- mcmapply(ki_prop_theta, ts_pred, MoreArgs =list(start_vec), mc.cores = num_cores)),
gcFirst = TRUE)[3]
)
theta_res_df <- data.frame("time" = ts_pred,
"Counts" = counts_theta_vec,
"Ki67high proportions" = ki_prop_theta_vec) %>%
gather(-time, key = "features", value = "results")
ggplot()+
geom_line(data = theta_res_df, aes(time, results), col = "#D87534", size = 1.5) +
labs(title = paste0("Source cohort - naive ", Population, " T cells"), x = "Host age (days)", y = NULL) +
facet_wrap( .~ features, scales = "free") + myTheme
print(time_integ_source)Cells that entered division
We assume that cells enter division with a rate \(\rho(a)\), depending on their age and then acquire the maximun ki67 intesity \(k = k_\text{max} = 1\), since \(k\) is the nomalised Ki67 intesity.
We also account for the delay-time that the cells that just entered division require before they renter the cell-cycle, by allowing cell division only in cells below a limiting \(k\) intensity \(k_\text{dL}\). The division limit i.e. \(k_\text{dL}\) lies between the thresold Ki67\(^{+}\) gating intensity and \(k_\text{max}\), such that \(k_\text{dL} \in (\bar{k}, 1)\).
Therefore the last boundary condition for the PDE in eq. , \[\begin{equation} u_\text{div}(t, a, k_\text{max} = 1) = \rho(a) \, u(t, a, k) dk, \quad \text{where } k \in (0, k_\text{dL}) \end{equation}\]
The system of ODEs for this charcteristic curve is, \[\begin{eqnarray} \begin{aligned} &\frac{dk}{ds} = -\beta \, k \quad \rightarrow k = k_0 \, e^{-\beta \, s} \quad \rightarrow k_0 =1 \quad \rightarrow s = \frac{-log(k)}{\beta} \\ &\frac{da}{ds} = 1 \quad \rightarrow a = s + a_0 \quad \rightarrow a_0 = a + \frac{log(k)}{\beta} \\ &\frac{dt}{ds} = 1 \quad \rightarrow t = s + t_0 \quad \rightarrow t_0 = t + \frac{log(k)}{\beta} \end{aligned} \end{eqnarray}\]
Using \(u_\text{div}(t_0, a_0, k_0), \, \rho(a)\) and \(\delta(a)\) the cohort of recently divided cells can be tracked on the characteristic curve `\(S\)‘, \[\begin{equation} u_\text{div}(t, a, k) = \rho(a) \, u(t + \frac{log(k)}{\beta}, a + \frac{log(k)}{\beta}, k) \, exp\big(- \int_{a + \frac{log(k)}{\beta}}^a \delta(\tau) \, d\tau \big) \end{equation}\] The distribution of ki67 expression within the cohort of cells that enetered division is, \[\begin{equation} u_\text{div}(t, k) = \int_0^t \rho(a) \, u(t + \frac{log(k)}{\beta}, a + \frac{log(k)}{\beta}, k) \, exp\big(- \int_{a + \frac{log(k)}{\beta}}^a \delta(\tau) \, d\tau \big) \, da. \end{equation}\] Therefore the total size (\(N_\text{d}\)) of the ’Division cohort’ and the porportion of Ki67\(^+\) cells (\(\kappa_\text{d}\)) within it are given by, \[\begin{eqnarray} \begin{aligned} & N_\text{div}(t) = \int_0^1 u_\text{d}(t, k) \; dk \\ & \kappa_\text{div}(t) = \frac{\int_{\bar{k}}^1 u_\text{d}(t, k) \, dk }{N_\text{d}(t)}. \end{aligned} \end{eqnarray}\]
## tracking cells that entered division
## age and ki67 distribution of the pool of cells that can enter division
N_pool_k_a <- function(age, ki, Time, parms){
## combined age and ki distribution of theta and init cohorts
value = N_theta_k_a(age, ki, Time, parms) + N_init_k_a(age, ki, Time, parms)
return(value)
}
### cellular division changing with cellage
rho_age <- function(age, parms){
rho0_log = parms[4] ### log transforming parameter ---> delta0 > 0
rho0 = exp(rho0_log)
r_rho = parms[5] ### free parameter
t0 = 5
return(rho0 * exp(-r_rho * age))
}
### recently divided cells k=1
N_rdiv_k_a <- function(age, ki, Time, parms){
delta0_log = parms[2] ### log transforming parameter ---> delta0 > 0
delta0 = exp(delta0_log)
r_del = parms[3] ### free parameter
t0 = 5
k_dl = 0.6
### pooled cells
u_pool_t0 <- function(ki, age , Time, parms){
rho0_log = parms[4] ### log transforming parameter ---> delta0 > 0
rho0 = exp(rho0_log)
rho0 * (N_init_a_k(ki, age + log(ki)/beta, Time + log(ki)/beta, parms) +
N_theta_a_k(ki, age + log(ki)/beta, Time + log(ki)/beta, parms))
}
N_rdiv_t0 <- function(age, ki, Time, parms){
delta0_log = parms[2] ### log transforming parameter ---> delta0 > 0
delta0 = exp(delta0_log)
r_del = parms[3] ### free parameter
t0 = 5
k_dl = 0.6
ifelse(ki <= k_dl,
u_pool_t0(ki, age, Time, parms) *
exp(- (delta0/r_del) * (exp(-r_del * (age + log(ki)/beta)) - exp(-r_del * age))),
0)
}
u_pool_t <- function(ki, age , Time, parms){
rho0_log = parms[4] ### log transforming parameter ---> delta0 > 0
rho0 = exp(rho0_log)
rho0 * (N_init_a_k(ki, age + log(ki)/beta, Time + log(ki)/beta, parms) +
N_theta_a_k(ki, age + log(ki)/beta, Time + log(ki)/beta, parms)+
N_rdiv_t0(age, ki, Time, parms))
}
N_rdiv_t <- function(age, ki, Time, parms){
delta0_log = parms[2] ### log transforming parameter ---> delta0 > 0
delta0 = exp(delta0_log)
r_del = parms[3] ### free parameter
t0 = 5
k_dl = 0.6
ifelse(ki <= k_dl,
u_pool_t(ki, age, Time, parms) *
exp(- (delta0/r_del) * (exp(-r_del * (age + log(ki)/beta)) - exp(-r_del * age))),
0)
# u_div_t(age, Time, parms) *
#exp(- (delta0/r_del) * (exp(-r_del * (age + log(ki)/beta)) - exp(-r_del * age)))
}
func_value <- function(age, ki, Time, Parms){
if(Time < t0 - (log(k_dl)/beta)){
0
} else if(Time == t0 - (log(k_dl)/beta)){
N_rdiv_t0(age, ki, Time, parms)
} else if (Time > t0 - (log(k_dl)/beta)){
N_rdiv_t(age, ki, Time,parms)
}
}
## mapply return the out put of the value function across each elements of vectors
## k and t for the same values of params as a list
return(mapply(func_value, age, ki, Time, MoreArgs = list(parms)))
}
## ki67 distribution of the cohorts that recently divided
N_rdiv_k <- function(ki, Time, parms){
rdiv_value <- function(ki, Time, parms){
t0 =5
return(integrate(N_rdiv_k_a, lower = 0, upper = Time - t0, ki=ki, Time=Time, parms = parms)$value)
}
## mapply return the out put of the value function across each elements of vectors
## k and t for the same values of params as a list
return(mapply(rdiv_value, ki, Time, MoreArgs = list(parms)))
}
## number of ki67+ cells within the cohorts that recently divided -- at time 't'
k_rdiv <- function(Time, parms){
k_value <- function(Time, parms){
k_bar = 1/exp(1)
## inetgrating across k_bar to 1
return(integrate(N_rdiv_k, lower = k_bar, upper = 1, Time=Time, parms=parms)$value)
}
## mapply return the out put of the value function across each elements of vector
## t for the same values of params as a list
return(mapply(k_value, Time, MoreArgs = list(parms)))
}
## total population size of the cohorts that recently divided -- at time 't'
N_rdiv <- function(Time, parms){
rdiv_value <- function(Time, parms){
## inetgrating across all k values
return(integrate(N_rdiv_k, lower = 0, upper = 1, Time=Time, parms=parms)$value)
}
## mapply return the out put of the value function across each elements of vector
##t for the same values of params as a list
return(mapply(rdiv_value, Time, MoreArgs = list(parms)))
}
## proportions of ki67+ cells within the cohorts that recently divided -- at time 't'
ki_prop_rdiv <- function(Time, parms){
# numbe rof ki67+ cells/ totla population size
value =
k_rdiv(Time, parms)/
N_rdiv(Time, parms)
return(value)
}# integration
time_rdiv_integ <- paste0("time taken for integration -- ",
system.time(
c(counts_rdiv_vec <- mcmapply(N_rdiv, ts_pred, MoreArgs =list(start_vec), mc.cores = num_cores),
ki_prop_rdiv_vec <- mcmapply(ki_prop_rdiv, ts_pred, MoreArgs =list(start_vec), mc.cores = num_cores)),
gcFirst = TRUE)[3]
)
rdiv_res_df <- data.frame("time" = ts_pred,
"Counts" = counts_rdiv_vec,
"Ki67high proportions" = ki_prop_rdiv_vec) %>%
gather(-time, key = "features", value = "results")
ggplot()+
geom_line(data = rdiv_res_df, aes(time, results), col = "#D87534", size = 1.5) +
labs(title = paste0("Source cohort - naive ", Population, " T cells"), x = "Host age (days)", y = NULL) +
facet_wrap( .~ features, scales = "free") + myTheme
print(time_rdiv_integ)The whole naive T cell compartment
Now, we can combine theese three cohorts to track the dynamics of the naive CD4 T cells pool size, which is an aggregate of the cells that have enetered from the thymus since \(t_0\), cells that were present at \(t_0\) and cells that were born from division. Also, the total poroportion of ki67\(^+\) cells is derived by summing the counts of ki67\(^+\) cells from all three cohorts and dividing it by total pool size.
\[\begin{eqnarray} \begin{aligned} &N(t) = N_\text{s}(t) + N_\text{p}(t) + N_\text{d}(t) \\ &\kappa(t) = \frac{\int_{\bar{k}}^1 u_\text{init}(t, k) + u_\text{source}(t, k) + u_\text{div}(t, k) \, dk }{N(t)}. \end{aligned} \end{eqnarray}\]
time_total_integ <- paste0("time taken for integration -- ",
system.time(
c(counts_total_vec <- mcmapply(N_total, ts_pred, MoreArgs =list(start_vec), mc.cores = num_cores),
ki_prop_total_vec <- mcmapply(ki_prop_total, ts_pred, MoreArgs =list(start_vec), mc.cores = num_cores)),
gcFirst = TRUE)[3])
total_res_df <- data.frame("time" = ts_pred,
"Counts" = counts_total_vec,
"Ki67high proportions" = ki_prop_total_vec) %>%
gather(-time, key = "features", value = "results")
ggplot()+
geom_line(data = total_res_df, aes(time, results), col = "#D87534", size = 1.5) +
labs(title = paste0("Source cohort - naive ", Population, " T cells"), x = "Host age (days)", y = NULL) +
facet_wrap( .~ features, scales = "free") + myTheme
print(time_total_integ)overlap_res_df <- data.frame("time" = ts_pred,
"Counts_theta" = counts_theta_vec,
"Counts_init" = counts_init_vec ,
"Counts_rdiv" = counts_rdiv_vec,
"Ki67high proportions" = ki_prop_total_vec) %>%
gather(-time, key = "features", value = "results")
ggplot()+
geom_line(data = total_res_df, aes(time, results), size = 1.5) +
labs(title = "Counts", x = "Time post t0 (days)", y = NULL) +
facet_wrap( .~ features, scales = "free")