Accepted Manuscript
Mathematical model of susceptibility, resistance, and resilience in the within-host dynamics between a Plasmodium parasite and the immune system Yi Yan, Brian Adam, Mary Galinski, Jessica Kissinger, Alberto Moreno, Juan B. Gutierrez PII: DOI: Reference:
S0025-5564(15)00209-6 10.1016/j.mbs.2015.10.003 MBS 7698
To appear in:
Mathematical Biosciences
Received date: Revised date: Accepted date:
26 February 2015 2 October 2015 7 October 2015
Please cite this article as: Yi Yan, Brian Adam, Mary Galinski, Jessica Kissinger, Alberto Moreno, Juan B. Gutierrez, Mathematical model of susceptibility, resistance, and resilience in the within-host dynamics between a Plasmodium parasite and the immune system, Mathematical Biosciences (2015), doi: 10.1016/j.mbs.2015.10.003
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights • We study the interaction between a Plasmodium parasite and the immune system.
CR IP T
• We propose a mathematical model for this interaction and prove that the model is well posed.
• We determine the conditions for resilience, resistance and susceptibility of the host.
AC
CE
PT
ED
M
AN US
• We conduct numerical simulations demonstrating the mains results.
1
ACCEPTED MANUSCRIPT
CR IP T
Mathematical model of susceptibility, resistance, and resilience in the within-host dynamics between a Plasmodium parasite and the immune system Yi Yan1 , Brian Adam2 , Mary Galinski4,5 , Jessica Kissinger1,3 , Alberto Moreno4 , Juan B. Gutierrez1,2,6
AN US
Abstract
PT
ED
M
We developed a coupled age-structured partial differential equation model to capture the disease dynamics during blood-stage malaria. The addition of age structure for the parasite population, with respect to previous models, allows us to better characterize the interaction between the malaria parasite and red blood cells during infection. Here we prove that the system we propose is well-posed and there exist at least two global states. We further demonstrate that the numerical simulation of the system coincides with clinically observed outcomes of primary and secondary malaria infection. The well-posedness of this system guarantees that the behavior of the model remains smooth, bounded, and continuously dependent on initial conditions; calibration with clinical data will constrain domains of parameters and variables to physiological ranges. 1. Introduction
AC
CE
The goal of this article is to formulate a mathematical model to characterize the within-host dynamics present in the disease malaria, between different host cell types and pathogens of the Plasmodium species. Plasmodium is a large genus of parasitic protozoa (unicellular eukaryotic organisms), 1
Institute of Bioinformatics, University of Georgia Department of Mathematics, University of Georgia 3 Department of Genetics, University of Georgia 4 Department of Medicine, Emory University 5 Department of Microbiology and Immunology, Emory University 6 Corresponding Author: [email protected], University of Georgia, Athens, GA 30602 2
Preprint submitted to Mathematical Biosciences
October 22, 2015
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
with complex genomes and sophisticated life cycles. The genome, behavior, and epidemiological characterization of Plasmodium are orders of magnitude more complex than that of viruses or bacteria. The clinical manifestation of infection by Plasmodium, the disease malaria, has a wide spectrum of symptoms, varying from asymptomatic to highly severe. Malaria affects birds, reptiles, and some mammals (mostly rodents and primates). The Plasmodium life cycle is comprised of several stages. The infection process in humans starts with the injection of sporozoites by female Anopheline mosquitoes into the skin of the host. This is followed by the liver stage, in which the inoculated sporozoites grow and multiply asexually within hepatocytes for 1-2 weeks to produce merozoites. The newly produced merozoites emerge from the liver and enter the blood stream. The blood-stage infection starts immediately after the hepatic stage; merozoites invade red blood cells (RBCs) where they also reproduce asexually. In some instances, the parasites develop into the sexual stage form called gametocytes. Gametocyteinfected red blood cells (iRBCs) infect newly feeding Anopheline mosquitoes and through sexual reproduction followed by many more rounds of asexual multiplication are ultimately transformed into sporozoites, thus completing the infection cycle between an Anopheline mosquito and its host. The parasite’s blood-stage infection in both human and non-human primates generally has a regular cycle of 24 to 72 hours depending on the species of the Plasmodium parasite [1, 2]. The parasites invade healthy RBCs and replicate asexually, remodeling and ultimately destroying the RBCs in the process. The destruction of RBCs during blood-stage malaria infection sometimes results in severe anemia, which is one of the major complications of malaria and a leading cause of mortality. Human and non-human primate RBCs have a normal life span of 120 to 100 days, respectively, where afterwards RBCs are cleared rapidly [3]. Both the age structure of RBCs and iRBCs play an important role in the hematodynamics of the host during malaria infection. It has been shown previously that Plasmodium vivax infects RBCs of different age groups at different rates and the interaction between the life cycle of the malaria parasite and the immune system can lead to synchrony of the parasite [4, 5]. Anemia caused by malaria infection is a complex process not only involving the malaria parasites and host RBCs, but also the innate and adaptive immune system of the host [6, 7]. It has been previously postulated that there are at least three factors that contribute to severe anemia resulting from blood stage malarial infection: (i) The destruction of RBCs directly by the parasite, (ii) the destruction of RBCs by a 3
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
non-parasitic factor, possibly the host immune system, and (iii) the disruption of the host’s erythropoietic process [6, 8, 3]; erythropoiesis refers to the differentiation and maturation of hematopoietic stem cells (HSCs) through committed cell lineages culminating in the production of erythrocytes/RBCs. Classic models of erythropoiesis rest on three variables representing the precursor/progenitor population of immature erythrocytes, the population of mature RBCs, and the concentration of erythropoietin (Epo), a signaling molecule critical to erythrocyte production. Mackey and Milton (1990) represented these three variables with a system of ordinary differential equations (ODEs) with one delay accounting for the time for Epo dependent maturation of progenitors into erythrocytes (about 6 days) [9]. In 1995, Belair et al. developed an “age-structured” partial differential equation (PDE) model for erythropoiesis [10]. Belair et al.’s model characterizes the proliferation and aging of the precursor and mature erythrocyte populations as nonlinear, first order PDEs and includes an ODE to represent Epo dynamics over time. Building on Belair’s 1995 model, Mahaffy et al. (1998) investigated the addition of state-dependent delays to the model [11]. Due to the difficulty involved in measuring Epo accurately over an extended period of time, we formulated the rate at which mature RBCs enter circulation as only dependent upon the current number of RBCs and the number of RBCs at equilibrium. Such a formulation has previously been shown by Savil et al. to capture the erythropoiesis process using data derived from a murine model [12]. As the only formulation of an erythropoiesis model without an Epo component that has been validated by time course data, we believe that such a formulation would also be able to capture the fluctuation in the erythropoietic process due to loss of RBCs within our hematodynamic system. Most recent and most relevant to our study, Thibodeaux et al. published two papers (2010, 2013) modifying Belair, Mackey, and Mahaffy’s 1995 model in an attempt to simulate and analyze erythropoietic dynamics subject to malaria infection [13, 14]. They conjoined an ODE-based model of malaria infection with an age-structured PDE model of erythropoiesis, and examined the dynamics of this new system to simulate the effects of hemozoin (Hz) on the suppression of erythropoiesis. In order to capture the complex interactions between different age groups of iRBCs, RBCs and various immune cells during blood stage malaria [15], we propose a coupled age-structured PDE model of both iRBCs and host RBCs. Our model borrows from the erythropoiesis PDE model [13, 14] and 4
ACCEPTED MANUSCRIPT
CR IP T
expands upon it by adding the parasite age structure system and the effect of immune cells on both RBCs and iRBCs. This paper is organized as follows: Section 2 contains the detailed derivation and description of our hematodynamic model. Section 3 proves the well posedness of the system and the existence of at least two global behaviors of biological interest. Section 4 provides numerical simulation results of our model. Section 5 offers some conclusions, and presents future directions. 2. Model Formulation
AC
CE
PT
ED
M
AN US
Let u(a, t) be a function which approximates the concentration of RBCs of age a at time point t. RBCs of age > amax are rapidly cleared from the circulation, thus we assume that u(a, t) = 0 when a > amax . Additionally, Let v(α, t) be a function that approximates the concentration of iRBCs of age α at time t. iRBCs of age αmax burst to produce merozoites, which in turn infect other RBCs, thus we assume that v(α, t) = 0 when α > αmax . Assuming that the loss of RBCs during blood stage malarial infection is only due to parasite invasion, innate immune cells, adaptive immune cells, random loss of RBCs during aging and the rapid clearance of RBCs with age > amax , and the destruction of RBCs by innate and adaptive immune cells follows the law of mass action, then the difference between the concentrations of RBCs across all ages between two time points, t1 and t2 is Z amax Z t2 Z t2 Z amax u(a, t2 ) da = u(a, t1 ) da + u(a0 , t) dt − u(amax , t) dt a0 t1 t1 a0 Z amax Z t2 v(αmax , t)r(a)p(u(a, t)) dt da −γ
where
−
− −
t1 a0 Z T X amax Z t2 i=1 Q
a0
XZ
wi (t)θi u(a, t) dt da
t1
amax Z t2
si (t)ψi u(a, t) dt da
t1 i=1 a0 Z amax Z t2
h(a)u(a, t) dt da ,
a0
(1)
t1
u(a, t) p(u(a, t)) = R amax u(a, t) a0 5
(2)
ACCEPTED MANUSCRIPT
t1
a0
a0
CR IP T
Ra is a probability density function such that a0max p(u(a, t))da = 1. The left hand side (LHS) term of (1) and the first three terms of the right hand side (RHS) of (1), Z t2 Z t2 Z amax Z amax u(amax , t) dt, u(a0 , t) dt − u(a, t1 ) da + u(a, t2 ) da = t1
a0
t1
AN US
(3) account for the change of the concentration of RBCs from t1 to t2 due to the production of new RBCs and the rapid clearance of RBCs with age > amax . The term Z amax Z t2 −γ v(αmax , t)r(a)p(u(a, t)) dt da , (4)
M
accounts for the loss of RBCs due to infection from bursting iRBCs. v(αmax , t) is the concentration of bursting iRBCs at time t. r(a) is the success rate of merozoites infecting RBCs of age a. γ is the expected number of merozoites each bursting iRBC produces. p(u(a, t)) is the expected percentage of merozoites infecting RBCs of specific age a at time t. The concentration of new iRBCs produced at time t is Z amax γ v(αmax , t)r(a)p(u(a, t)) da. (5) The term
ED
a0
PT
−
T Z X i=1
amax
a0
Z
t2
wi (t)θi u(a, t) dt da,
(6)
t1
AC
CE
accounts for the loss of RBCs due to T different kinds of innate immune cells. wi (t) is the concentration of the ith kind of innate immune cells at time t. θi represents the ith kind of innate immune cell’s effectiveness at destroying RBCs. The term Q Z amax Z t2 X si (t)ψi u(a, t) dt da, (7) i=1
a0
t1
accounts for the loss of RBCs due to Q different kinds of adaptive immune cells. si (t) is the concentration of the ith kind of adaptive immune cells at time t and ψi represents the ith kind of adaptive immune cell’s effectiveness at destroying RBCs. 6
ACCEPTED MANUSCRIPT
The last term
amax Z t2
Z
h(a)u(a, t) dt da ,
(8)
t1
a0
Z
t2
t1
u(a0 , t) dt −
Z
t2
u(amax , t) dt = −
t1
Z
t1
t2 Z amax
ua da dt.
AN US
and
a0
a0
a0
CR IP T
accounts for the random loss of RBCs during aging. h(a) is the natural death rate of RBCs of different ages. Using the fundamental theorem of Calculus, it is clear that: Z amax Z t2 Z amax Z amax ut dt da, (9) u(a, t1 ) da = u(a, t2 ) da −
t1
(10)
a0
Substituting (9) and (10) into (1), we obtain: Z amax Z t2 Z t2 Z amax ut dt da + ua da dt = a0 t1 t1 a0 Z amax Z t2 −γ v(αmax , t)r(a)p(u(a, t)) dt da i=1
amax Z t2
a0
amax Z t2
a0
i=1 amax Z t2
Z
PT
−
Q Z X
−
a0
wi (t)θi u(a, t) dt da
t1
ED
−
t1
M
a0
T Z X
si (t)ψi u(a, t) dt da
t1
h(a)u(a, t) dt da,
(11)
t1
AC
CE
which can be rearrange to the following equation because all the terms have the same integral: Z amax Z t2 (ut + ua ) da dt = a0 t1 Z amax Z t2 −γ v(αmax , t)r(a)p(u(a, t))dt da a0 t1 ! Z amax Z t2 X Q T X − wi (t)θi + si (t)ψi + h(a) u(a, t)dt da. (12) a0
t1
i=1
i=1
7
ACCEPTED MANUSCRIPT
− γv(αmax , t)r(a)p(u(a, t)).
(13) is subject to the boundary condition:
AN US
u(0, t) = f (t, ϕ(t)),
CR IP T
If the function u(a, t) and its partial derivatives are continuous functions then the integrands on both sides of (12) should be equal, which leads us to the following partial differential equation: ! Q T X X ∂u ∂u si (t)ψi + h(a) u(a, t) + =− wi (t)θi + ∂t ∂a i=1 i=1
where
( ς, t < Td f (t, ϕ(t)) = ςeε(ϕ0 −ϕ(t)) , t > Td .
(13)
(14)
(15)
AC
CE
PT
ED
M
This boundary condition reflects the erythropoietic response to the change in concentration of RBCs. ϕ(t) is the concentration of RBCs at time point t, ς is the normal rate at which RBCs enter the peripheral blood, Td reflects a one time lag in the erythropoietic response to the change in concentration of RBCs, as the erythropoietic process only responds to the infection after a time period of Td . ϕ0 is the normal concentration of RBCs in a healthy host and ε is a parameter that controls the maximum amount of erythropoietic response. During malarial infection, as the concentration of RBCs decreases due to malaria parasites, the rate of new RBCs entering into circulation increases. The age structured model of iRBCs can be similarly derived. Let v(α, t) be a function that approximates the concentration of iRBCs of age α at time t. iRBCs of age αmax bursts to produce merozoites, which in turn infect other RBCs, thus we assume that v(α, t) = 0 when α > αmax . Assuming the loss of iRBCs is only due to the bursting of iRBCs of age αmax , innate immune cells and adaptive immune cells, then the difference between the concentration of
8
ACCEPTED MANUSCRIPT
iRBCs across all ages between two time points, t1 and t2 is Z Z t2 Z αmax Z αmax v(α0 , t) dt − v(α, t1 ) dα + v(α, t2 ) dα = − −
T Z X
αmax Z t2
i=1 Q
α0
i=1
α0
XZ
wi (t)φi v(α, t) dt dα
t1
αmax Z t2
v(αmax , t) dt
t1
t1
CR IP T
α0
α0
t2
si (t)bi (t)v(α, t) dt dα
t1
(16)
AN US
Using the same technique we used to derive (13), we arrive at the following PDE for v(α, t): ! Q T X X ∂v ∂v +V =− wi (t)φi + si (t)bi (t) v(α, t), (17) ∂t ∂α i=1 i=1
ED
M
where V is the speed at which the parasite ages, which is subject to the unit of α and the term νi , (18) bi (t) = Rt 1 + exp −λi ( t0 v(0, t) dt − ξi )
AC
CE
PT
describes the adaptive immune cell i’s effector strength against iRBCs. The function bi (t) captures the ability of adaptive immune cells to increase their effect against iRBCs through increased exposure, which increases, as the Rt exposure to new iRBCs, characterized by t0 v(0, t) increases past threshold ξi . The effector strength has a maximum value of νi as the amount of exposure to iRBCs becomes substantially large. The term ! T X − wi (t)φi v(α, t), (19) i=1
describes the destruction of iRBC by different innate immune cells each with constant strength φi . The term ! Q X − si (t)bi (t) v(α, t), (20) i=1
9
ACCEPTED MANUSCRIPT
CR IP T
describes the destruction of iRBCs by different adaptive immune cells each with strength bi (t) that changes according to the total exposure to iRBCs. Additionally, this PDE is subject to the boundary condition: Z amax v(αmax , t)r(a)p(u(a, t))da, (21) v(α0 , t) = γ a0
and oi (t) = $i +
AN US
which describes the creation of new iRBCs. Because the life span of an iRBC is between 24 and 72 hours, and the infection process is within minutes, we assume that such process is instantaneous within our model. Additionally, due to the presence of iRBCs, the following equation is proposed to describe the change in concentration of innate immune cells: Z αmax dwi = oi (t) − βi wi (t) − τi wi v(α, t)dα, (22) dt α0 i − $i . R αmax 1 + exp −ηi ( α0 v(α, t)dα − Mi )
(23)
CE
PT
ED
M
We assume that each innate immune cell population w(t)i decays at a rate of βi . The production of innate immune cells increases as the concentration of iRBCs increases past a threshold Mi and reaches a maximum level of i − $i . Additionally, the term τi describes the loss of functionality of innate immune cells upon contact with iRBCs. (23) describes the change in the rate of production of innate immune cells in response to the change in concentration of iRBCs. The change in adaptive immune cell population in response to increased concentration of iRBCs is similarly formulated as the following: Z αmax dsi = li (t) − δi si (t) − ϑi si v(α, t)dα, (24) dt α0
AC
and
li (t) = σi +
%i − σi . R αmax 1 + exp −ωi ( α0 v(α, t)dα − Ri )
(25)
We assume that each adaptive immune cell population s(t)i decays at a rate of δi . The production of adaptive immune cells increases as the concentration of iRBCs increases past a threshold Ri and reaches a maximum level %i − 10
ACCEPTED MANUSCRIPT
− γv(αmax , t)r(a)p(u(a, t)), T X
wi (t)φi +
i=1
Q X
!
si (t)bi (t) v(α, t),
AN US
∂v ∂v +V =− ∂t ∂α
CR IP T
σi . Additionally, the term ϑi describes the loss of functionality of adaptive immune cells upon contact with iRBCs. (25) describes the change in the the rate of production of adaptive immune cells in response to the change in concentration of iRBCs. Putting (13), (14), (17), (21), (22), and (24) together we arrive at the following PDE system: ! Q T X X ∂u ∂u si (t)ψi + h(a) u(a, t) + =− wi (t)θi + ∂t ∂a i=1 i=1
i=1
(26)
Z αmax dwi =oi (t) − βi wi (t) − τi wi (t) v(α, t)dα, dt α0 Z αmax dsi =li (t) − δi si (t) − ϑi si (t) v(α, t)dα, dt α0
M
subject to the following initial and boundary conditions:
PT
ED
u(a, 0) =g(a), u(0, t) =f (t, ϕ(t)), v(α, 0) =c(α), Z amax v(0, t) =γ v(αmax , t)r(a)p(u(a, t)) da,
(27)
a0
AC
CE
where g(a) is the initial RBC age distribution and c(α) is the initial iRBC age distribution. Since only Plasmodium vivax has been shown to infect RBCs of different ages at different rates [5], we consider r(a) to be a constant κ within our model where κ ∈ [0, 1]. Furthermore, due to a lack of data regarding the loss of innate and adaptive immune function upon interaction with iRBC, we consider ϑi and τi to be 0 within our model. Additionally, according to the erythropoiesis model developed by Savil et al [12], h(a) was also set to be a constant. Lastly, we assume that the total number of iRBCs is much smaller than the total number of RBCs. Therefore, the crowding effect on the production of iRBCs is negligible. This assumption is justified if the host 11
ACCEPTED MANUSCRIPT
− γκv(αmax , t)p(u(a, t)),
∂v ∂v +V =− ∂t ∂α
T X
wi (t)φi +
i=1
Q X i=1
!
si (t)bi (t) v(α, t),
AN US
dwi =oi (t) − βi wi (t), dt dsi =li (t) − δi si (t), dt
CR IP T
becomes resilient or completely removes the infection. The simplified model is: ! Q T X X ∂u ∂u si (t)ψi + h(a) u(a, t) + =− wi (t)θi + ∂t ∂a i=1 i=1
(28)
subject to the following initial and boundary conditions:
M
u(a, 0) =g(a), u(0, t) =f (t, ϕ(t)), v(α, 0) =c(α), Z amax v(0, t) =γκ v(αmax , t)p(u(a, t)) da.
(29)
ED
a0
PT
From this point on we will be working with the above simplified version of our model. 3. System Behavior
CE
3.1. Immune Saturation Given an large parasite exposure, the adaptive immune strength bi (t) becomes constant:
AC
Rt
t0
lim
v(0,t) dt→+∞
νi Rt v(0,t) dt→+∞ 1 + exp −λi ( t0 v(0, t) dt − ξi ) t0
bi (t) = R t
lim
(30)
=νi .
Additionally, the production rate of both innate and adaptive immune cells also reaches a constant value given a large iRBC population in the 12
ACCEPTED MANUSCRIPT
system: v(a,t)dα→+∞
α0
oi (t) = R αmax lim
v(a,t)dα→+∞
α0
+
=i ,
$i
i − $i , R αmax 1 + exp −ηi ( α0 v(α, t)dα − Mi )
and α0
lim
v(a,t)dα→+∞
li (t) = R αmax lim
v(a,t)dα→+∞
α0
+
=%i .
σi
%i − σi , R αmax 1 + exp −ωi ( α0 v(α, t)dα − Ri )
AN US
R αmax
(31)
CR IP T
lim
R αmax
(32)
M
Because the production rate of both the innate and adaptive immune cells becomes constant given a large parasite population, the population of both the innate and adaptive immune cells also become constant: R αmax
v(a,t)dα→+∞
ED
α0
lim
PT
and
R αmax α0
lim
v(a,t)dα→+∞
wi (t) =
i βi
(33)
si (t) =
%i . δi
(34)
AC
CE
Assuming immune saturation by substituting (33) and (34),(13) becomes: ! Q T X X i %i ∂u ∂u + =− θi + ψi + h(a) u(a, t) ∂t ∂a βi δi (35) i=1 i=1 − γκv(αmax , t)p(u(a, t)),
where h(a) is a constant decay term describing the random loss of RBCs over time and i , βi , θi , %i , δi , ψi are all constants. The system further simplifies to: ∂u ∂u + = C2 u(a, t) − γκv(αmax , t)p(u(a, t)). (36) ∂t ∂a
13
ACCEPTED MANUSCRIPT
CR IP T
Additionally, assuming immune saturation by substituting (33),(34) and (30), (16) becomes ! Q T X X %i ∂v ∂v i +V =− φi + νi v(α, t), (37) ∂t ∂α βi δi i=1 i=1
AN US
where i , βi , φi , %i , δi and νi are all positive constants. (37) further simplifies to: ∂v ∂v +V = C1 v(α, t). (38) ∂t ∂α Corollary 3.1. There exist an unique solution for v(α, t) under immune saturation conditions. Proof. Using the fact that α is a function of time, v(α, t) can be expressed as v(α(t), t), which transforms (38) into the following Equation: d ∂v dα ∂v v= + = C1 v, dt ∂t dt ∂α
ED
M
which reduces (38) to the following ODE system: ( dα = V, dt dv = C1 v, dt
(40) dv dt
= C1 v. Solving the
(41)
CE
PT
which states that along the curves given by dα = V, dt above system, we arrive at: ˆ (t) = V t + a0 , α vˆ(t) = v(tˆ0 ) exp(C1 t), vˆ(t0 ) = v(a0 , 0).
(39)
AC
However, due to the fact that our system is bounded by t ∈ [0, ∞] and α ˆ ∈ [0, αmax ], a0 does not exist when t > α(t) . But the boundary condition of V v(0, t) allows us to express (41) as the following ( t0 = t − α(t) , V (42) vˆ(t) = v(0, t0 ) exp(C1 (t − t0 )), which essentially traces the characteristic curve to the boundary v(0, t0 ) if t > α(t) . Consider that the solution to v(α, t) is a system of ODEs with V 14
ACCEPTED MANUSCRIPT
ln 1/γκ αmax /V of αmax . V
Corollary 3.2. If C1 = behavior with a period C1 >
ln 1/γκ αmax /V
CR IP T
= −Cv, which is Lipschitz continuous. Thus using the Picardthe form dv dt Lindel¨of theorem, each ODEs within the system have an unique solution, therefore, there exist an unique solution to the PDE v(α, t). , then the boundary v(0, t) exhibits periodic If C1
1 αmax v(0, t − ( )), γκ V 15
(49)
ACCEPTED MANUSCRIPT
and if C1 >
v(αmax , t)
0. Which leads to: sup u(0, t) ≤ ςeεϕ0 .
(54)
t
inf u(0, t) ≥ ςeε(ϕ0 −amax ςe t
εϕ0 )
>0
(55)
AC
Therefore u(0, t) ∈ (0, ςeεϕ0 ]. Corollary 3.5. There exists a set of parameters and initial conditions such that 0 < inf a,t u(a, t) < inf t u(0, t).
16
ACCEPTED MANUSCRIPT
inf uˆ(0) + C2
Z
amax
uˆ(t)dt − γκ
0
where p(ˆ u, t) =
ˆ(t) R amaxu u(a,t)da 0
inf C2 1 , amax
t0 +amax
t0
v(αmax , t)p(ˆ u, t)dt ≥ N, (56)
and C2 < 0. Assume u(a, t) ≥ N > 0:
Z
amax
0
Let |C2 |
0.
inf (1 + C2 amax )ˆ u(0) − γκ
Z
t0 +amax
t0
v(αmax , t)p(ˆ u, t)dt ≥ N
(58)
Z
− γκ
t0 +amax
t0
v(αmax , t)p(ˆ u, t)dt ≥ N − (1 + C2 amax ) inf uˆ(0) (59)
ED
inf
M
Let uˆ(0) = inf uˆ(0) and 0 < N < (1 + C2 amax ) inf uˆ(0) < inf uˆ(0), then:
Because both sides of the inequality are negative, we can rearrange the inequality in the following form: Z
v(αmax , t)p(ˆ u, t)dt ≤ (1 + C2 amax ) inf uˆ(0) − N
t0 +amax
PT
inf γκ
t0
(60)
AC
CE
ˆ(0) tu Because supt p(ˆ u, t) < sup = M and let v(αmax , t) either be oscillating or amax N decreasing as shown in Corollary 3.2, then
sup γκ
Z
t0 +amax
t0
v(αmax , t)p(ˆ u, t)dt ≤ γκamax
sup t0 ≤t≤t0 +amax
v(αmax , t)M (61)
Let supt0 ≤t≤t0 +amax v(αmax , t) = Vmax , then
sup γκ
Z
t0 +amax
t0
v(αmax , t)p(ˆ u, t)dt ≤ γκamax Vmax M. 17
(62)
ACCEPTED MANUSCRIPT
Divide (60) by γκamax M: (1 + C2 amax ) inf uˆ(0) − N γκamax M
Therefore, if Vmax ≤
(63)
CR IP T
Vmax ≤
(1 + C2 amax ) inf uˆ(0) − N , γκamax M 1 |C2 | < amax
then
0 < N < inf (u(a, t)) < inf (u(0, t)) .
t
AN US
a,t
Corollary 3.6. There exist an unique solution of u(a, t).
Proof. As stated in Corollary 3.3 the solution of u(a, t) is a system of ODEs with the form:
M
dˆ u = C2 uˆ − γκv(αmax , t)p(ˆ u, t) = F (ˆ u, t). dt
(64)
Let uˆ1 and uˆ2 be different, then:
CE
PT
ED
|F (ˆ u1 , t) − F (ˆ u2 , t)| = |C2 (ˆ u1 − uˆ2 ) − γκv(αmax , t)(p(ˆ u1 , t) − p(ˆ u2 , t))| . (65) R amax u ˆn Because p(ˆ un , t) = En (t) , where En (t) = a0 un (a, t)da. (65) further simplifies to: uˆ uˆ2 1 − . (66) |F (ˆ u1 , t) − F (ˆ u2 , t)| = C2 (ˆ u1 − uˆ2 ) − γκv(αmax , t) E1 (t) E2 (t)
AC
As a result, it follows that: u ˆ E (t) − u ˆ E (t) 1 2 2 1 |F (ˆ u1 , t) − F (ˆ u2 , t)| = C2 (ˆ u1 − uˆ2 ) − γκv(αmax , t) E1 (t)E2 (t) uˆ1 E2 (t) − uˆ2 E1 (t) ≤ |C2 ||ˆ u1 − uˆ2 | + γκv(αmax , t) E1 (t)E2 (t) uˆ E (t) − uˆ E (t) 1 2 2 1 ≤ |C2 ||ˆ u1 − uˆ2 | + |γκv(αmax , t)| E1 (t)E2 (t) (67) 18
ACCEPTED MANUSCRIPT
CR IP T
In the previous Corollary, we have proved that under suitable initial conditions and parameters, 0 < N < inf t u(a, t) < u(a, t) < supt u(a, t). Which implies that En (t) is bounded. Let A = inf t En (t) and B = supt En (t). Without loss of generality. assume E1 (t) ≤ E2 (t). It follows that (67) is bounded above by: supt E2 (t) |ˆ u1 − uˆ2 | inf t E1 (t) inf t E2 (t) t B ≤ |C2 ||ˆ u1 − uˆ2 | + sup |γκv(αmax , t)| 2 |ˆ u1 − uˆ2 | A B u1 − uˆ2 |. = |C2 | + sup |γκv(αmax , t)| 2 |ˆ A
AN US
sup |C2 ||ˆ u1 − uˆ2 | + sup |γκv(αmax , t)|
If we define:
B Ω = |C2 | + sup |γκv(αmax , t)| 2 , A then we have established that:
(69)
(70)
M
|F (ˆ u1 , t) − F (ˆ u2 , t)| ≤ Ω|ˆ u1 − uˆ2 |.
(68)
ED
We have demonstrated that each of the infinite number of ODEs composing the solution of u(a, t) is Lipschitz continuous. Applying the Picard-Lindel¨of theorem, we arrive at the conclusion that there exist an unique solution for each of the ODEs composing the solution of u(a, t) and therefore there exist an unique solution of u(a, t).
PT
Corollary 3.7. The unique solution for v(α, t) is stable for a set of parameters and initial conditions.
AC
CE
Proof. Assume that there are two solutions of the system v1 (α, t) and v2 (α, t) with different boundary conditions:v1 (0, t),v1 (α, 0),v2 (0, t) and v2 (α, 0). Recall Corollary 3.1 which states that the solution of v(α, t) is a system of ODEs with the following form: ˆ (t) = V t + a0 α (71) vˆ(t) = v(tˆ0 ) exp(C1 t) vˆ(t0 ) = v(a0 , 0). Using (41), it is clear that, v1 (α, t) = v1 (t − V t, 0) exp(C1 t) on its characteristic curves and v2 (α, t) = v2 (t − V t, 0) exp(C1 t) on its characteristic curves. 19
ACCEPTED MANUSCRIPT
Therefore, the difference between v1 (α, t) and v2 (α, t) at any given α and t can be expressed as the following if t < αmax : (72)
CR IP T
|v1 (α, t) − v2 (α, t)| = |(v1 (t − V t, 0) − v2 (t − V t, 0)) exp(C1 t)| .
If t > αmax , using (42), the difference between v1 (α, t) and v2 (α, t) at any given α and t can be expressed as the following: |v1 (α, t) − v2 (α, t)| = |(v1 (0, t0 ) − v2 (0, t0 )) exp(C1 (t − t0))| , where
(73)
α . (74) V Recall Corollary 3.2, which states that there exist a set of parameters such that v(α, t) is bounded. Therefore, according to (72) and (73) it is clear that the solution of v(α, t) depends continuously on its initial and boundary conditions, thus v(α, t) is stable.
AN US
t0 = t −
Corollary 3.8. The unique solution for u(a, t) is stable.
AC
CE
PT
ED
M
Proof. Assume that there are two solutions of the system u1 (a, t) and u2 (a, t) with different boundary conditions: u1 (0, t), u1 (a, 0), u2 (0, t) and u2 (a, 0). Then the difference between the total numbers of RBCs in u1 (a, t) and u2 (a, t) at an arbitrary time point t can be expressed in the following form as demon-
20
ACCEPTED MANUSCRIPT
amax
u1 (a, t1 ) − u2 (a, t1 ) da
a0
a0
+
Z
t2
t Z 1t2
CR IP T
strated in (1): Z Z amax u1 (a, t2 ) − u2 (a, t2 ) da =
u1 (a0 , t) − u2 (a0 , t) dt
u1 (amax , t) − u2 (amax , t) dt Z amax Z t2 v(αmax , t)r(a)p(u1 (a, t)) dt da −γ t1 a0 Z amax Z t2 +γ v(αmax , t)r(a)p(u2 (a, t)) dt da −
a0 t1 Z t2 Z T a max X i=1 Q
t1
amax Z t2
a0
i=1 amax Z t2
Z
a0
t1
t1
wi (t)θi u1 (a, t) − u2 (a, t) dt da
si (t)ψi u(a, t) − u2 (a, t) dt da
h(a)u1 (a, t) − u2 (a, t) dt da .
ED
−
a0
XZ
M
−
t1
AN US
−
AC
CE
PT
Let W = u1 (a, t) − u2 (a, t), then (75) simplifies to: Z amax Z amax Z t2 W (a, t2 ) da = W (a, t1 ) da + W (a0 , t) dt a0 a0 t1 Z t2 − W (amax , t) dt t1 Z amax Z t2 −γ v(αmax , t)r(a)p(u1 (a, t)) dt da a0 t1 Z amax Z t2 +γ v(αmax , t)r(a)p(u2 (a, t)) dt da −
a0 t1 Z Z t2 T a max X i=1
a0
t1
21
wi (t)θi W (a, t) dt da
(75)
ACCEPTED MANUSCRIPT
−
amax Z t2
si (t)ψi W (a, t) dt da
t1 i=1 a0 Z amax Z t2
h(a)W (a, t) dt da .
t1
a0
Assuming r(a) is a constant,then the following is true: Z t2 Z amax Z t2 v(αmax , t)r(a)dt. v(αmax , t)r(a)p(u1 (a, t)) dt da = γ γ a0
t1
t1
(76)
CR IP T
−
Q Z X
(77)
t1 T X Z amax Z t2 i=1 Q
a0
XZ
M
−
AN US
Substituting (77) into (76) we arrive at the following system: Z amax Z amax Z t2 W (a, t2 ) da = W (a, t1 ) da + W (a0 , t) dt a0 a0 t1 Z t2 − W (amax , t) dt
ED
−
−
wi (t)θi W (a, t) dt da
t1
amax Z t2
si (t)ψi W (a, t) dt da
t1 i=1 a0 Z amax Z t2
h(a)W (a, t) dt da .
a0
(78)
t1
CE
PT
If W (a, t) and its partial derivatives are continuous functions then (78) can be rearranged into an transport system as demonstrated in (36): ∂W ∂W + = C2 W (a, t), ∂t ∂a
(79)
AC
subject to the following initial and boundary condition: W (a, 0) =u1 (a, 0) − u2 (a, 0), W (0, t) =u1 (0, t) − u2 (0, t).
(80)
As shown in Corollary 3.1,the solution to (79) is a system of ODEs with the following form: ( da = 1, dt (81) ˆ dW ˆ, = C2 W dt
22
ACCEPTED MANUSCRIPT
which states that alone the characteristic curve given by The solution to (81) along the characteristic curve is:
da dW , dt dt
= exp(C2 t).
ˆ (t) = W ˆ (0) exp(C2 t), W
(82)
CR IP T
ˆ (0) is the initial condition for the ODE on the characteristic curve where W and lies either on W (a, 0) or W (0, a). Recall Corollary 3.5 and Corollary 3.4, which state that there exists a set of parameters such that u(a, t) is bounded. Therefore, according to (81) and (82), the solution to u(a, t) depends continuously on its initial and boundary conditions, thus u(a, t) is stable.
AN US
Lemma 3.9. The PDE system (26) is well-posed for a set of parameters and initial conditions.
M
Proof. Corollary 3.1 and Corollary 3.3 demonstrate that the solution of the two transport equation within (26) are systems of first order ODEs. Corollary 3.6 and Corollary 3.1 proved the existence and uniqueness of the solution to the first order ODEs composing the solution to the two transport equation within our system. Additionally, Corollary 3.7 and Corollary 3.8 demonstrate that the solution to (26) depends continuously on the initial and boundary conditions. Therefore, the PDE system (26) is well-posed.
ED
4. Numerical Solution
CE
where
PT
Recall that our system is a transport system with the form:
m(a, t) =
T X i=1
∂u ∂u + = −m(a, t)u(a, t), ∂t ∂a
wi (t)θi +
Q X i=1
γv(αmax , t)r(a) si (t)ψi + h(a) + R amax u(a, t) a0
(83) !
.
(84)
AC
Using the first order forward finite difference scheme at u(a, t), where the time interval is divided into N steps of size k (indicated with the superindex j), and the age interval is divided into M segments of size h (indicated with the subindex i). it is clear that: ∂u uj+1 − uji = i , ∂t k 23
(85)
ACCEPTED MANUSCRIPT
and
AN US
CR IP T
uj − uji−1 ∂u = i . (86) ∂a h Substituting the above Equation into (83), and let r = k/h, the system becomes: uji − uji−1 uj+1 − uji i + = −mji uji , (87) k h resulting in the following difference equation: j+1 j j i j ui = −kmj ui + (1 − r)ui + rui−1 , (88) u0i = g(i), j u0 = f (j). Using the above schema one can iteratively approximate all u(a, t) at t = k then t = 2k etc.
Lemma 4.1. j The above finite difference schema is numerically stable if r < 1 and max{ mi : ∀i, j}k < (2 − 2r).
PT
ED
M
Proof. Let uˆji be the numerical solution of (88). Then the error eji = uˆji − uji satisfies: j+1 j j i j ei = −kmj ei + (1 − r)ei + rei−1 , (89) e0i = 0, j e0 = 0. Let ej = max{ eji : ∀i} for each j ≥ 0. Then
AC
CE
j+1 e = (1 − r + kmij ) + rej i i−1 (90) i ≤ ( 1 − r + kmj + r)ej Thus ej+1 ≤ ( 1 − r − kmij + r)ej for all j ≥ 0. Let m ˆ = max{ mji : ∀i, j} for each i ≥ 0 and j ≥ 0. For r < 1 and k m ˆ < (2 − 2r), we have ej+1 ≤ ej for all j ≥ 0. This means that the error ej is not increasing when j is increasing. Therefore, (88) is stable. The same schema is also applied to v(α, t). Recall that: ∂v ∂v +V = −n(α, t)v(α, t), ∂t ∂α 24
(91)
ACCEPTED MANUSCRIPT
where n(α, t) =
T X
wi (t)φi +
Q X i=1
i=1
!
si (t)bi (t) .
(92)
CR IP T
Using the same finite difference schema one can iteratively approximate all v(α, t) at t = k then t = 2k etc. j+1 j j i j vi = −knj vi + (1 − rV )vi + rV vi−1 (93) vi0 = c(i) j v0 = d(j)
AN US
Which is also stable for rV < 1 and max{nij }k < (2 − 2rV ).
M
Additionally, using Taylor expansion, the value of w(t) and s(t) are approximated as following: dwi wi (t + k) = w(t) + k dt (t), R αmax dwi (94) v(a, t)dα, (t) = o (t) − β w (t) − τ w i i i i i dt α0 wi (0) = A,
PT
5. Results
dsi si (t + k) = s(t) + k dt (t), Rα dsi (t) = li (t) − δi si (t) − ϑi si α0max v(a, t)dα, dt si (0) = B.
(95)
ED
and
AC
CE
5.1. RBC Equilibrium Condition The erythropoietic process responds to events altering the concentration of RBCs and ensures that the number and age distribution of RBCs return to normal level upon perturbation. Figure 1 shows that using an arbitrary initial age distribution of RBCs, the system eventually returns to the basal RBC age distribution overtime. This result demonstrates that our system reflects the actual biological property of erythropoietic control, where the host responds to a perturbation in concentration of RBCs by altering accordingly, the rate of RBC production in the bone marrow, and the rate of RBCs entering circulation. This regulatory effect can be seen with the fluctuation of u(0, t) in Figure 1, which first responds to the lower than normal number of RBCs with a increase of RBCs entering circulation after a delay of 20 days. 25
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Figure 1: Using an arbitrary initial condition, the system returns to Equilibrium condition after 80 days.
AC
CE
PT
ED
5.2. Susceptibility As stated in Corollary 3.2, if the immune strength is weak, eventually the iRBC population grows without bounds and RBCs are depleted, resulting in severe anemia (Figure 2). Depletion of RBCs is characteristic of first time malaria infection which without treatment could result in mortality of the host. Within our model, this behavior can be achieved under at least two parameter sets. The first one is where C1 , the combined effector strength upon saturation of both innate and adaptive immune cells is weak, which allows for the depletion of RBCs. The second set is where C1 is large enough to halt the exponential growth of iRBCs, but the parameters controlling the speed at which innate immune cells, adaptive immune cells and adaptive immune effector strength growth are small, resulting in the depletion of RBCs before C1 can reach its maximum value. The second scenario reflects the case where a host is susceptible to severe anemia during primary infection, but during repeated infections, C1 had enough time to reach maximum value, and the host is resilient to severe anemia [3].
26
,
CR IP T
ACCEPTED MANUSCRIPT
AN US
Figure 2: The iRBC population grows exponentially due to weak immune strength
ED
M
5.3. Resilience
ln 1/γκ αmax /V
PT
Figure 3: iRBC and RBC co-exist when total immune strength is close to
AC
CE
1/γκ , the iRBC and RBC If the immune effector strength is close to αlnmax /V population coexist for a long period of time as shown in Figure 3. This state of resilience is characteristic of repeated malaria infection, where the host has already built up an immune strength against the parasite. However, sustained periodic coexistence is almost never observed in the real world, which is expected due to the fixed maximum immune strength required for sustained coexistence.
5.4. Resistance As stated in Corollary 3.2, if the immune strength reaches a greater value 1/γκ than αlnmax , then iRBC population decays constantly until the complete /V 27
CR IP T
ACCEPTED MANUSCRIPT
ln 1/γκ αmax /V
AN US
Figure 4: iRBC and RBC co-exist when total immune strength is greater than
removal of iRBC, as seen in Figure 4. This state of resistance is rarely observed without medical intervention during first infection; however, it is feasible to develop an immune response against malaria upon continuous and repeated exposure high enough to clear the parasite [16].
M
6. Discussion
AC
CE
PT
ED
In this paper we present a theoretical model to capture the complex interaction between iRBCs, RBCs and the host immune system during the course of a Plasmodium infection. Despite not knowing many of the parameters regarding the immune system’s effect on RBCs and iRBCs, our model was able to capture the general behavior of two major outcomes of infection with Plasmodium: (i) the depletion of host RBCs in a naive individual and (ii) the coexistence of RBCs, iRBCs, and the immune system in individuals with previous exposure. Our model states that the immune effector has to be equal to an exact number for sustained coexistence of the system, which on first look seems rather impossible as biological systems are rarely precise. This further coincides with known malaria infection dynamics where the oscillation of iRBC populations is rarely sustained or regular. To produce sustained oscillation over weeks or months, the immune effector strength only needs to lie in a range rather than being fixed. It is worth mentioning, that one of the major simplifications of the model was the assumption of r(a), the infection probability of RBCs of different age, to be constant. For sufficiently regular r(a), the well-posedness results can be derived in a similar fashion. However, the rigid trichotomy which 28
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
conveniently characterizes the system is not inherited. The simplified system in (28) results in uncoupled dynamics for iRBCs. This could imply that infected cells may appear even in the absence of uninfected cells; however, this condition would not be physiological. Given realistic initial and boundary conditions, the simplified system permits analysis and approximates well the reality of a malaria infection. We plan to calibrate our model in the future using time course data from the Malaria Host-Pathogen Interaction Center (MaHPIC), which is a large collaborative project between Emory University, The University of Georgia, The Georgia Institute of Technology, and the Centers for Disease Control and Prevention. The MaHPIC project aims to obtain time courses of -omic data (transcriptomic, metabolomic, lipidomic, proteomic) from non-human primates infected with several different species of Plasmodium parasites during multiple experiments, with initial experiments lasting for about 100 days. Using data produced by the MaHPIC group, which include daily RBC counts, iRBC counts and immune cell counts over the 100-day period, we plan to calibrate the full model, including the different infection rates of different age groups of RBCs and the different effector strengths against different age groups of iRBCs. We expect to find discrepancies between our model prediction and the actual experimental data; such discrepancy will not only provide biological insight regarding the hematodynamics of malaria infection, but will also allow us to re-iterate our model. Our current PDE model is by no means a complete reflection of the complex within-host dynamics of malaria due to its vast simplification of the immune system’s effect on the iRBC and RBC populations along with its complete disregard of the within immune system interaction. However, the model could be expanded to also contain these features and serve as a liaison where the molecular features of the system can be reflected at the cellular level. We plan to use the molecular level data generated by MaHPIC to further expand our model to incorporate the more subtle aspects of the within-host interaction of malaria infection. Acknowledgments The authors would like to thank Jacob A. Grey for his help in refining the corollaries 3.4 through 3.7, and also Stacey Lapp, Tracey Lamb, Rabindra Tirouvanziam, Chet Joyner, and Esmeralda Meyer for fruitful discussions. Two anonymous reviewers helped to improve this manuscript significantly. 29
ACCEPTED MANUSCRIPT
CR IP T
This research has been supported Malaria Host-Pathogen Interaction Center - MaHPIC, NIH’s NIAID contract HHSN272201200031C, and the International Centers for Excellence in Malaria Research - Center for non-Amazonian regions of Latin America (CLAIM), NIH’s NIAID cooperative agreement U19AI089702. References
[1] Z. Bozdech, M. Llin´as, B. L. Pulliam, E. D. Wong, J. Zhu, J. L. DeRisi, The transcriptome of the intraerythrocytic developmental cycle of plasmodium falciparum, PLoS biology 1 (1) (2003) e5.
AN US
[2] W. E. Collins, M. Warren, J. S. Sullivan, G. G. Galland, Plasmodium coatneyi: observations on periodicity, mosquito infection, and transmission to macaca mulatta monkeys., The American journal of tropical medicine and hygiene 64 (3) (2001) 101–110.
M
[3] A. Moreno, M. Cabrera-Mora, A. Garcia, J. Orkin, E. Strobert, J. W. Barnwell, M. R. Galinski, Plasmodium coatneyi in rhesus macaques replicates the multisystemic dysfunction of severe malaria in humans, Infection and immunity 81 (6) (2013) 1889–1904.
ED
[4] I. M. Rouzine, F. E. McKenzie, Link between immune response and parasite synchronization in malaria, Proceedings of the National Academy of Sciences 100 (6) (2003) 3473–3478.
PT
[5] D. A. Moreno-P´erez, J. A. Ru´ız, M. A. Patarroyo, Reticulocytes: Plasmodium vivax target cells, Biology of the Cell 105 (6) (2013) 251–260.
CE
[6] D. J. Perkins, T. Were, G. C. Davenport, P. Kempaiah, J. B. Hittner, J. M. Ong’echa, Severe malarial anemia: innate immunity and pathogenesis, International journal of biological sciences 7 (9) (2011) 1427.
AC
[7] T. J. Lamb, J. Langhorne, The severity of malarial anaemia in plasmodium chabaudi infections of balb/c mice is determined independently of the number of circulating parasites, Malaria journal 7 (1) (2008) 68. [8] N. M. Douglas, N. M. Anstey, P. A. Buffet, J. R. Poespoprodjo, T. W. Yeo, N. J. White, R. N. Price, et al., The anaemia of plasmodium vivax malaria, Malar J 11 (135.10) (2012) 1186. 30
ACCEPTED MANUSCRIPT
[9] M. Mackey, J. Milton, Feedback, delays and the origin of blood cell dynamics, University of Minnesota. Institute for Mathematics and Its Applications, 1990.
CR IP T
[10] J. Belair, M. C. Mackey, J. M. Mahaffy, Age-structured and two-delay models for erythropoiesis, Mathematical Biosciences 128 (1) (1995) 317– 346. [11] J. M. Mahaffy, J. Belair, M. C. Mackey, Hematopoietic model with moving boundary condition and state dependent delay] applications in erythropoiesis, Journal of Theoretical Biology 190 (2) (1998) 135–146.
AN US
[12] N. J. Savill, W. Chadwick, S. E. Reece, Quantitative analysis of mechanisms that govern red blood cell age structure and dynamics during anaemia, PLoS computational biology 5 (6) (2009) e1000416. [13] J. J. Thibodeaux, Modeling erythropoiesis subject to malaria infection, Mathematical biosciences 225 (1) (2010) 59–67.
M
[14] A. S. Ackleh, B. Ma, J. J. Thibodeaux, A second-order high resolution finite difference scheme for a structured erythropoiesis model subject to malaria infection, Mathematical biosciences.
PT
ED
[15] C. Metcalf, A. Graham, S. Huijben, V. Barclay, G. Long, B. Grenfell, A. Read, O. Bjørnstad, Partitioning regulatory mechanisms of withinhost malaria dynamics using the effective propagation number, Science 333 (6045) (2011) 984–988.
AC
CE
[16] D. L. Doolan, C. Dobano, J. K. Baird, Acquired immunity to malaria, Clinical microbiology reviews 22 (1) (2009) 13–36.
31
ACCEPTED MANUSCRIPT
APENDIX: Table of Parameters and Variables Description
si (t)
h(a) f (t, ϕt ) r(a)
AC
CE
PT
ED
li (t)
M
p(u(a, t))
oi (t)
Source
n/a n/a n/a
n/a n/a n/a
n/a n/a n/a n/a
n/a n/a n/a n/a
n/a
n/a
0.02 0.03
MaHPIC
AN US
u(a, t) v(α, t) ϕ(t) wi (t)
bi (t)
Value
CR IP T
t a α
Unit Independent Variables Time. Day Age of RBC. Day Age of iRBC. 10− 1Day Dependent Variables Circulating RBC age distribution. cells/ul/day iRBC age distribution. cells/ul/day Circulating RBC concentration. cells/ul Circulating innate immune cell concencells/ul tration. Circulating adaptive immune cell con- cells/ul centration. Functions Percentage of RBC that leaves circula1/day tion per day, usually a constant. Rate at which RBC enters circulation. cells/ul/day Success rate of merozoite invading 1/day RBC of age. a Probability of a merozoite infecting a dimensionless RBC of age. a at time t Rate at which adaptive immunity cell 1/(day · cells/ul) i destroys iRBC. Rate at which innate immune cell encells/ul/day ters circulation. Rate at which adaptive immune cell cells/ul/day enters circulation.
32
n/a [0, 1]
n/a Literature
n/a
n/a
n/a
Estimated
n/a
Estimated
n/a
Estimated
ACCEPTED MANUSCRIPT
ψi φi βi δi τi ϑi λi $i σi i %i νi ηi ωi Mi
CE
Ri
xii
AC
ε ς
Td
Value
Source
10
n/a
6400000
MaHPIC
25
Literature
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
CR IP T
θi
M
γ
ED
ϕ0
PT
V
Unit Parameters Speed at which parasite ages, equal to 1/(day) 1/α. Normal circulating RBC concentracells/ul tion, equal to ϕ(0). Average number of merozoite pro- dimensionless duced by a single iRBC. Rate at which innate immunity cell i 1/(day · cells/ul) destroys RBC. Rate at which adaptive immunity cell 1/(day · cells/ul) i destroys RBC. Rate at which innate immunity cell i 1/(day · cells/ul) destroys iRBC. Rate at which innate immunity cell de- 1/day cays. Rate at which adaptive immunity cell 1/day decays. Rate at which innate immunity cell de- 1/(day · cells/ul) cay due to contact with iRBC. Rate at which adaptive immunity cell 1/(day · cells/ul) decay due to contact with iRBC. Coefficient for change of adaptive im- 1/cells/ul munity effector strength. Normal rate of production of innate cells/ul/day immune cell i. Normal rate of production of adaptive cells/ul/day immune cell i. Maximum rate of production of innate cells/ul/day immune cell i. Maximum rate of production of adapcells/ul/day tive immune cell i. Maximum adaptive immunity effector 1/day/cells/ul i strength. Coefficient for change of production of 1/cells/ul innate immune cell i. Coefficient for change of production of 1/cells/ul adaptive immune cell i. Amount of parasite where the produc- cells/ul tion of innate immune cell i. increase the most Amount of parasite where the produc- cells/ul tion of adaptive immune cell i. increase the most Amount of parasite where the produc- cells/ul tion of adaptive immune cell i. effector strength increase the most Coefficient for change of erythro- 1/cells/ul poiesis. Number of RBC entering blood cells/ul/day stream. Delay of hematopoietic response. day
AN US
Description
33
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
8000
MaHPIC
10 20
Literature
Mathematical model of susceptibility, resistance, and resilience in the within-host dynamics between a Plasmodium parasite and the immune system Yi Yan, Brian Adam, Mary Galinski, Jessica Kissinger, Alberto Moreno, Juan B. Gutierrez PII: DOI: Reference:
S0025-5564(15)00209-6 10.1016/j.mbs.2015.10.003 MBS 7698
To appear in:
Mathematical Biosciences
Received date: Revised date: Accepted date:
26 February 2015 2 October 2015 7 October 2015
Please cite this article as: Yi Yan, Brian Adam, Mary Galinski, Jessica Kissinger, Alberto Moreno, Juan B. Gutierrez, Mathematical model of susceptibility, resistance, and resilience in the within-host dynamics between a Plasmodium parasite and the immune system, Mathematical Biosciences (2015), doi: 10.1016/j.mbs.2015.10.003
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights • We study the interaction between a Plasmodium parasite and the immune system.
CR IP T
• We propose a mathematical model for this interaction and prove that the model is well posed.
• We determine the conditions for resilience, resistance and susceptibility of the host.
AC
CE
PT
ED
M
AN US
• We conduct numerical simulations demonstrating the mains results.
1
ACCEPTED MANUSCRIPT
CR IP T
Mathematical model of susceptibility, resistance, and resilience in the within-host dynamics between a Plasmodium parasite and the immune system Yi Yan1 , Brian Adam2 , Mary Galinski4,5 , Jessica Kissinger1,3 , Alberto Moreno4 , Juan B. Gutierrez1,2,6
AN US
Abstract
PT
ED
M
We developed a coupled age-structured partial differential equation model to capture the disease dynamics during blood-stage malaria. The addition of age structure for the parasite population, with respect to previous models, allows us to better characterize the interaction between the malaria parasite and red blood cells during infection. Here we prove that the system we propose is well-posed and there exist at least two global states. We further demonstrate that the numerical simulation of the system coincides with clinically observed outcomes of primary and secondary malaria infection. The well-posedness of this system guarantees that the behavior of the model remains smooth, bounded, and continuously dependent on initial conditions; calibration with clinical data will constrain domains of parameters and variables to physiological ranges. 1. Introduction
AC
CE
The goal of this article is to formulate a mathematical model to characterize the within-host dynamics present in the disease malaria, between different host cell types and pathogens of the Plasmodium species. Plasmodium is a large genus of parasitic protozoa (unicellular eukaryotic organisms), 1
Institute of Bioinformatics, University of Georgia Department of Mathematics, University of Georgia 3 Department of Genetics, University of Georgia 4 Department of Medicine, Emory University 5 Department of Microbiology and Immunology, Emory University 6 Corresponding Author: [email protected], University of Georgia, Athens, GA 30602 2
Preprint submitted to Mathematical Biosciences
October 22, 2015
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
with complex genomes and sophisticated life cycles. The genome, behavior, and epidemiological characterization of Plasmodium are orders of magnitude more complex than that of viruses or bacteria. The clinical manifestation of infection by Plasmodium, the disease malaria, has a wide spectrum of symptoms, varying from asymptomatic to highly severe. Malaria affects birds, reptiles, and some mammals (mostly rodents and primates). The Plasmodium life cycle is comprised of several stages. The infection process in humans starts with the injection of sporozoites by female Anopheline mosquitoes into the skin of the host. This is followed by the liver stage, in which the inoculated sporozoites grow and multiply asexually within hepatocytes for 1-2 weeks to produce merozoites. The newly produced merozoites emerge from the liver and enter the blood stream. The blood-stage infection starts immediately after the hepatic stage; merozoites invade red blood cells (RBCs) where they also reproduce asexually. In some instances, the parasites develop into the sexual stage form called gametocytes. Gametocyteinfected red blood cells (iRBCs) infect newly feeding Anopheline mosquitoes and through sexual reproduction followed by many more rounds of asexual multiplication are ultimately transformed into sporozoites, thus completing the infection cycle between an Anopheline mosquito and its host. The parasite’s blood-stage infection in both human and non-human primates generally has a regular cycle of 24 to 72 hours depending on the species of the Plasmodium parasite [1, 2]. The parasites invade healthy RBCs and replicate asexually, remodeling and ultimately destroying the RBCs in the process. The destruction of RBCs during blood-stage malaria infection sometimes results in severe anemia, which is one of the major complications of malaria and a leading cause of mortality. Human and non-human primate RBCs have a normal life span of 120 to 100 days, respectively, where afterwards RBCs are cleared rapidly [3]. Both the age structure of RBCs and iRBCs play an important role in the hematodynamics of the host during malaria infection. It has been shown previously that Plasmodium vivax infects RBCs of different age groups at different rates and the interaction between the life cycle of the malaria parasite and the immune system can lead to synchrony of the parasite [4, 5]. Anemia caused by malaria infection is a complex process not only involving the malaria parasites and host RBCs, but also the innate and adaptive immune system of the host [6, 7]. It has been previously postulated that there are at least three factors that contribute to severe anemia resulting from blood stage malarial infection: (i) The destruction of RBCs directly by the parasite, (ii) the destruction of RBCs by a 3
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
non-parasitic factor, possibly the host immune system, and (iii) the disruption of the host’s erythropoietic process [6, 8, 3]; erythropoiesis refers to the differentiation and maturation of hematopoietic stem cells (HSCs) through committed cell lineages culminating in the production of erythrocytes/RBCs. Classic models of erythropoiesis rest on three variables representing the precursor/progenitor population of immature erythrocytes, the population of mature RBCs, and the concentration of erythropoietin (Epo), a signaling molecule critical to erythrocyte production. Mackey and Milton (1990) represented these three variables with a system of ordinary differential equations (ODEs) with one delay accounting for the time for Epo dependent maturation of progenitors into erythrocytes (about 6 days) [9]. In 1995, Belair et al. developed an “age-structured” partial differential equation (PDE) model for erythropoiesis [10]. Belair et al.’s model characterizes the proliferation and aging of the precursor and mature erythrocyte populations as nonlinear, first order PDEs and includes an ODE to represent Epo dynamics over time. Building on Belair’s 1995 model, Mahaffy et al. (1998) investigated the addition of state-dependent delays to the model [11]. Due to the difficulty involved in measuring Epo accurately over an extended period of time, we formulated the rate at which mature RBCs enter circulation as only dependent upon the current number of RBCs and the number of RBCs at equilibrium. Such a formulation has previously been shown by Savil et al. to capture the erythropoiesis process using data derived from a murine model [12]. As the only formulation of an erythropoiesis model without an Epo component that has been validated by time course data, we believe that such a formulation would also be able to capture the fluctuation in the erythropoietic process due to loss of RBCs within our hematodynamic system. Most recent and most relevant to our study, Thibodeaux et al. published two papers (2010, 2013) modifying Belair, Mackey, and Mahaffy’s 1995 model in an attempt to simulate and analyze erythropoietic dynamics subject to malaria infection [13, 14]. They conjoined an ODE-based model of malaria infection with an age-structured PDE model of erythropoiesis, and examined the dynamics of this new system to simulate the effects of hemozoin (Hz) on the suppression of erythropoiesis. In order to capture the complex interactions between different age groups of iRBCs, RBCs and various immune cells during blood stage malaria [15], we propose a coupled age-structured PDE model of both iRBCs and host RBCs. Our model borrows from the erythropoiesis PDE model [13, 14] and 4
ACCEPTED MANUSCRIPT
CR IP T
expands upon it by adding the parasite age structure system and the effect of immune cells on both RBCs and iRBCs. This paper is organized as follows: Section 2 contains the detailed derivation and description of our hematodynamic model. Section 3 proves the well posedness of the system and the existence of at least two global behaviors of biological interest. Section 4 provides numerical simulation results of our model. Section 5 offers some conclusions, and presents future directions. 2. Model Formulation
AC
CE
PT
ED
M
AN US
Let u(a, t) be a function which approximates the concentration of RBCs of age a at time point t. RBCs of age > amax are rapidly cleared from the circulation, thus we assume that u(a, t) = 0 when a > amax . Additionally, Let v(α, t) be a function that approximates the concentration of iRBCs of age α at time t. iRBCs of age αmax burst to produce merozoites, which in turn infect other RBCs, thus we assume that v(α, t) = 0 when α > αmax . Assuming that the loss of RBCs during blood stage malarial infection is only due to parasite invasion, innate immune cells, adaptive immune cells, random loss of RBCs during aging and the rapid clearance of RBCs with age > amax , and the destruction of RBCs by innate and adaptive immune cells follows the law of mass action, then the difference between the concentrations of RBCs across all ages between two time points, t1 and t2 is Z amax Z t2 Z t2 Z amax u(a, t2 ) da = u(a, t1 ) da + u(a0 , t) dt − u(amax , t) dt a0 t1 t1 a0 Z amax Z t2 v(αmax , t)r(a)p(u(a, t)) dt da −γ
where
−
− −
t1 a0 Z T X amax Z t2 i=1 Q
a0
XZ
wi (t)θi u(a, t) dt da
t1
amax Z t2
si (t)ψi u(a, t) dt da
t1 i=1 a0 Z amax Z t2
h(a)u(a, t) dt da ,
a0
(1)
t1
u(a, t) p(u(a, t)) = R amax u(a, t) a0 5
(2)
ACCEPTED MANUSCRIPT
t1
a0
a0
CR IP T
Ra is a probability density function such that a0max p(u(a, t))da = 1. The left hand side (LHS) term of (1) and the first three terms of the right hand side (RHS) of (1), Z t2 Z t2 Z amax Z amax u(amax , t) dt, u(a0 , t) dt − u(a, t1 ) da + u(a, t2 ) da = t1
a0
t1
AN US
(3) account for the change of the concentration of RBCs from t1 to t2 due to the production of new RBCs and the rapid clearance of RBCs with age > amax . The term Z amax Z t2 −γ v(αmax , t)r(a)p(u(a, t)) dt da , (4)
M
accounts for the loss of RBCs due to infection from bursting iRBCs. v(αmax , t) is the concentration of bursting iRBCs at time t. r(a) is the success rate of merozoites infecting RBCs of age a. γ is the expected number of merozoites each bursting iRBC produces. p(u(a, t)) is the expected percentage of merozoites infecting RBCs of specific age a at time t. The concentration of new iRBCs produced at time t is Z amax γ v(αmax , t)r(a)p(u(a, t)) da. (5) The term
ED
a0
PT
−
T Z X i=1
amax
a0
Z
t2
wi (t)θi u(a, t) dt da,
(6)
t1
AC
CE
accounts for the loss of RBCs due to T different kinds of innate immune cells. wi (t) is the concentration of the ith kind of innate immune cells at time t. θi represents the ith kind of innate immune cell’s effectiveness at destroying RBCs. The term Q Z amax Z t2 X si (t)ψi u(a, t) dt da, (7) i=1
a0
t1
accounts for the loss of RBCs due to Q different kinds of adaptive immune cells. si (t) is the concentration of the ith kind of adaptive immune cells at time t and ψi represents the ith kind of adaptive immune cell’s effectiveness at destroying RBCs. 6
ACCEPTED MANUSCRIPT
The last term
amax Z t2
Z
h(a)u(a, t) dt da ,
(8)
t1
a0
Z
t2
t1
u(a0 , t) dt −
Z
t2
u(amax , t) dt = −
t1
Z
t1
t2 Z amax
ua da dt.
AN US
and
a0
a0
a0
CR IP T
accounts for the random loss of RBCs during aging. h(a) is the natural death rate of RBCs of different ages. Using the fundamental theorem of Calculus, it is clear that: Z amax Z t2 Z amax Z amax ut dt da, (9) u(a, t1 ) da = u(a, t2 ) da −
t1
(10)
a0
Substituting (9) and (10) into (1), we obtain: Z amax Z t2 Z t2 Z amax ut dt da + ua da dt = a0 t1 t1 a0 Z amax Z t2 −γ v(αmax , t)r(a)p(u(a, t)) dt da i=1
amax Z t2
a0
amax Z t2
a0
i=1 amax Z t2
Z
PT
−
Q Z X
−
a0
wi (t)θi u(a, t) dt da
t1
ED
−
t1
M
a0
T Z X
si (t)ψi u(a, t) dt da
t1
h(a)u(a, t) dt da,
(11)
t1
AC
CE
which can be rearrange to the following equation because all the terms have the same integral: Z amax Z t2 (ut + ua ) da dt = a0 t1 Z amax Z t2 −γ v(αmax , t)r(a)p(u(a, t))dt da a0 t1 ! Z amax Z t2 X Q T X − wi (t)θi + si (t)ψi + h(a) u(a, t)dt da. (12) a0
t1
i=1
i=1
7
ACCEPTED MANUSCRIPT
− γv(αmax , t)r(a)p(u(a, t)).
(13) is subject to the boundary condition:
AN US
u(0, t) = f (t, ϕ(t)),
CR IP T
If the function u(a, t) and its partial derivatives are continuous functions then the integrands on both sides of (12) should be equal, which leads us to the following partial differential equation: ! Q T X X ∂u ∂u si (t)ψi + h(a) u(a, t) + =− wi (t)θi + ∂t ∂a i=1 i=1
where
( ς, t < Td f (t, ϕ(t)) = ςeε(ϕ0 −ϕ(t)) , t > Td .
(13)
(14)
(15)
AC
CE
PT
ED
M
This boundary condition reflects the erythropoietic response to the change in concentration of RBCs. ϕ(t) is the concentration of RBCs at time point t, ς is the normal rate at which RBCs enter the peripheral blood, Td reflects a one time lag in the erythropoietic response to the change in concentration of RBCs, as the erythropoietic process only responds to the infection after a time period of Td . ϕ0 is the normal concentration of RBCs in a healthy host and ε is a parameter that controls the maximum amount of erythropoietic response. During malarial infection, as the concentration of RBCs decreases due to malaria parasites, the rate of new RBCs entering into circulation increases. The age structured model of iRBCs can be similarly derived. Let v(α, t) be a function that approximates the concentration of iRBCs of age α at time t. iRBCs of age αmax bursts to produce merozoites, which in turn infect other RBCs, thus we assume that v(α, t) = 0 when α > αmax . Assuming the loss of iRBCs is only due to the bursting of iRBCs of age αmax , innate immune cells and adaptive immune cells, then the difference between the concentration of
8
ACCEPTED MANUSCRIPT
iRBCs across all ages between two time points, t1 and t2 is Z Z t2 Z αmax Z αmax v(α0 , t) dt − v(α, t1 ) dα + v(α, t2 ) dα = − −
T Z X
αmax Z t2
i=1 Q
α0
i=1
α0
XZ
wi (t)φi v(α, t) dt dα
t1
αmax Z t2
v(αmax , t) dt
t1
t1
CR IP T
α0
α0
t2
si (t)bi (t)v(α, t) dt dα
t1
(16)
AN US
Using the same technique we used to derive (13), we arrive at the following PDE for v(α, t): ! Q T X X ∂v ∂v +V =− wi (t)φi + si (t)bi (t) v(α, t), (17) ∂t ∂α i=1 i=1
ED
M
where V is the speed at which the parasite ages, which is subject to the unit of α and the term νi , (18) bi (t) = Rt 1 + exp −λi ( t0 v(0, t) dt − ξi )
AC
CE
PT
describes the adaptive immune cell i’s effector strength against iRBCs. The function bi (t) captures the ability of adaptive immune cells to increase their effect against iRBCs through increased exposure, which increases, as the Rt exposure to new iRBCs, characterized by t0 v(0, t) increases past threshold ξi . The effector strength has a maximum value of νi as the amount of exposure to iRBCs becomes substantially large. The term ! T X − wi (t)φi v(α, t), (19) i=1
describes the destruction of iRBC by different innate immune cells each with constant strength φi . The term ! Q X − si (t)bi (t) v(α, t), (20) i=1
9
ACCEPTED MANUSCRIPT
CR IP T
describes the destruction of iRBCs by different adaptive immune cells each with strength bi (t) that changes according to the total exposure to iRBCs. Additionally, this PDE is subject to the boundary condition: Z amax v(αmax , t)r(a)p(u(a, t))da, (21) v(α0 , t) = γ a0
and oi (t) = $i +
AN US
which describes the creation of new iRBCs. Because the life span of an iRBC is between 24 and 72 hours, and the infection process is within minutes, we assume that such process is instantaneous within our model. Additionally, due to the presence of iRBCs, the following equation is proposed to describe the change in concentration of innate immune cells: Z αmax dwi = oi (t) − βi wi (t) − τi wi v(α, t)dα, (22) dt α0 i − $i . R αmax 1 + exp −ηi ( α0 v(α, t)dα − Mi )
(23)
CE
PT
ED
M
We assume that each innate immune cell population w(t)i decays at a rate of βi . The production of innate immune cells increases as the concentration of iRBCs increases past a threshold Mi and reaches a maximum level of i − $i . Additionally, the term τi describes the loss of functionality of innate immune cells upon contact with iRBCs. (23) describes the change in the rate of production of innate immune cells in response to the change in concentration of iRBCs. The change in adaptive immune cell population in response to increased concentration of iRBCs is similarly formulated as the following: Z αmax dsi = li (t) − δi si (t) − ϑi si v(α, t)dα, (24) dt α0
AC
and
li (t) = σi +
%i − σi . R αmax 1 + exp −ωi ( α0 v(α, t)dα − Ri )
(25)
We assume that each adaptive immune cell population s(t)i decays at a rate of δi . The production of adaptive immune cells increases as the concentration of iRBCs increases past a threshold Ri and reaches a maximum level %i − 10
ACCEPTED MANUSCRIPT
− γv(αmax , t)r(a)p(u(a, t)), T X
wi (t)φi +
i=1
Q X
!
si (t)bi (t) v(α, t),
AN US
∂v ∂v +V =− ∂t ∂α
CR IP T
σi . Additionally, the term ϑi describes the loss of functionality of adaptive immune cells upon contact with iRBCs. (25) describes the change in the the rate of production of adaptive immune cells in response to the change in concentration of iRBCs. Putting (13), (14), (17), (21), (22), and (24) together we arrive at the following PDE system: ! Q T X X ∂u ∂u si (t)ψi + h(a) u(a, t) + =− wi (t)θi + ∂t ∂a i=1 i=1
i=1
(26)
Z αmax dwi =oi (t) − βi wi (t) − τi wi (t) v(α, t)dα, dt α0 Z αmax dsi =li (t) − δi si (t) − ϑi si (t) v(α, t)dα, dt α0
M
subject to the following initial and boundary conditions:
PT
ED
u(a, 0) =g(a), u(0, t) =f (t, ϕ(t)), v(α, 0) =c(α), Z amax v(0, t) =γ v(αmax , t)r(a)p(u(a, t)) da,
(27)
a0
AC
CE
where g(a) is the initial RBC age distribution and c(α) is the initial iRBC age distribution. Since only Plasmodium vivax has been shown to infect RBCs of different ages at different rates [5], we consider r(a) to be a constant κ within our model where κ ∈ [0, 1]. Furthermore, due to a lack of data regarding the loss of innate and adaptive immune function upon interaction with iRBC, we consider ϑi and τi to be 0 within our model. Additionally, according to the erythropoiesis model developed by Savil et al [12], h(a) was also set to be a constant. Lastly, we assume that the total number of iRBCs is much smaller than the total number of RBCs. Therefore, the crowding effect on the production of iRBCs is negligible. This assumption is justified if the host 11
ACCEPTED MANUSCRIPT
− γκv(αmax , t)p(u(a, t)),
∂v ∂v +V =− ∂t ∂α
T X
wi (t)φi +
i=1
Q X i=1
!
si (t)bi (t) v(α, t),
AN US
dwi =oi (t) − βi wi (t), dt dsi =li (t) − δi si (t), dt
CR IP T
becomes resilient or completely removes the infection. The simplified model is: ! Q T X X ∂u ∂u si (t)ψi + h(a) u(a, t) + =− wi (t)θi + ∂t ∂a i=1 i=1
(28)
subject to the following initial and boundary conditions:
M
u(a, 0) =g(a), u(0, t) =f (t, ϕ(t)), v(α, 0) =c(α), Z amax v(0, t) =γκ v(αmax , t)p(u(a, t)) da.
(29)
ED
a0
PT
From this point on we will be working with the above simplified version of our model. 3. System Behavior
CE
3.1. Immune Saturation Given an large parasite exposure, the adaptive immune strength bi (t) becomes constant:
AC
Rt
t0
lim
v(0,t) dt→+∞
νi Rt v(0,t) dt→+∞ 1 + exp −λi ( t0 v(0, t) dt − ξi ) t0
bi (t) = R t
lim
(30)
=νi .
Additionally, the production rate of both innate and adaptive immune cells also reaches a constant value given a large iRBC population in the 12
ACCEPTED MANUSCRIPT
system: v(a,t)dα→+∞
α0
oi (t) = R αmax lim
v(a,t)dα→+∞
α0
+
=i ,
$i
i − $i , R αmax 1 + exp −ηi ( α0 v(α, t)dα − Mi )
and α0
lim
v(a,t)dα→+∞
li (t) = R αmax lim
v(a,t)dα→+∞
α0
+
=%i .
σi
%i − σi , R αmax 1 + exp −ωi ( α0 v(α, t)dα − Ri )
AN US
R αmax
(31)
CR IP T
lim
R αmax
(32)
M
Because the production rate of both the innate and adaptive immune cells becomes constant given a large parasite population, the population of both the innate and adaptive immune cells also become constant: R αmax
v(a,t)dα→+∞
ED
α0
lim
PT
and
R αmax α0
lim
v(a,t)dα→+∞
wi (t) =
i βi
(33)
si (t) =
%i . δi
(34)
AC
CE
Assuming immune saturation by substituting (33) and (34),(13) becomes: ! Q T X X i %i ∂u ∂u + =− θi + ψi + h(a) u(a, t) ∂t ∂a βi δi (35) i=1 i=1 − γκv(αmax , t)p(u(a, t)),
where h(a) is a constant decay term describing the random loss of RBCs over time and i , βi , θi , %i , δi , ψi are all constants. The system further simplifies to: ∂u ∂u + = C2 u(a, t) − γκv(αmax , t)p(u(a, t)). (36) ∂t ∂a
13
ACCEPTED MANUSCRIPT
CR IP T
Additionally, assuming immune saturation by substituting (33),(34) and (30), (16) becomes ! Q T X X %i ∂v ∂v i +V =− φi + νi v(α, t), (37) ∂t ∂α βi δi i=1 i=1
AN US
where i , βi , φi , %i , δi and νi are all positive constants. (37) further simplifies to: ∂v ∂v +V = C1 v(α, t). (38) ∂t ∂α Corollary 3.1. There exist an unique solution for v(α, t) under immune saturation conditions. Proof. Using the fact that α is a function of time, v(α, t) can be expressed as v(α(t), t), which transforms (38) into the following Equation: d ∂v dα ∂v v= + = C1 v, dt ∂t dt ∂α
ED
M
which reduces (38) to the following ODE system: ( dα = V, dt dv = C1 v, dt
(40) dv dt
= C1 v. Solving the
(41)
CE
PT
which states that along the curves given by dα = V, dt above system, we arrive at: ˆ (t) = V t + a0 , α vˆ(t) = v(tˆ0 ) exp(C1 t), vˆ(t0 ) = v(a0 , 0).
(39)
AC
However, due to the fact that our system is bounded by t ∈ [0, ∞] and α ˆ ∈ [0, αmax ], a0 does not exist when t > α(t) . But the boundary condition of V v(0, t) allows us to express (41) as the following ( t0 = t − α(t) , V (42) vˆ(t) = v(0, t0 ) exp(C1 (t − t0 )), which essentially traces the characteristic curve to the boundary v(0, t0 ) if t > α(t) . Consider that the solution to v(α, t) is a system of ODEs with V 14
ACCEPTED MANUSCRIPT
ln 1/γκ αmax /V of αmax . V
Corollary 3.2. If C1 = behavior with a period C1 >
ln 1/γκ αmax /V
CR IP T
= −Cv, which is Lipschitz continuous. Thus using the Picardthe form dv dt Lindel¨of theorem, each ODEs within the system have an unique solution, therefore, there exist an unique solution to the PDE v(α, t). , then the boundary v(0, t) exhibits periodic If C1
1 αmax v(0, t − ( )), γκ V 15
(49)
ACCEPTED MANUSCRIPT
and if C1 >
v(αmax , t)
0. Which leads to: sup u(0, t) ≤ ςeεϕ0 .
(54)
t
inf u(0, t) ≥ ςeε(ϕ0 −amax ςe t
εϕ0 )
>0
(55)
AC
Therefore u(0, t) ∈ (0, ςeεϕ0 ]. Corollary 3.5. There exists a set of parameters and initial conditions such that 0 < inf a,t u(a, t) < inf t u(0, t).
16
ACCEPTED MANUSCRIPT
inf uˆ(0) + C2
Z
amax
uˆ(t)dt − γκ
0
where p(ˆ u, t) =
ˆ(t) R amaxu u(a,t)da 0
inf C2 1 , amax
t0 +amax
t0
v(αmax , t)p(ˆ u, t)dt ≥ N, (56)
and C2 < 0. Assume u(a, t) ≥ N > 0:
Z
amax
0
Let |C2 |
0.
inf (1 + C2 amax )ˆ u(0) − γκ
Z
t0 +amax
t0
v(αmax , t)p(ˆ u, t)dt ≥ N
(58)
Z
− γκ
t0 +amax
t0
v(αmax , t)p(ˆ u, t)dt ≥ N − (1 + C2 amax ) inf uˆ(0) (59)
ED
inf
M
Let uˆ(0) = inf uˆ(0) and 0 < N < (1 + C2 amax ) inf uˆ(0) < inf uˆ(0), then:
Because both sides of the inequality are negative, we can rearrange the inequality in the following form: Z
v(αmax , t)p(ˆ u, t)dt ≤ (1 + C2 amax ) inf uˆ(0) − N
t0 +amax
PT
inf γκ
t0
(60)
AC
CE
ˆ(0) tu Because supt p(ˆ u, t) < sup = M and let v(αmax , t) either be oscillating or amax N decreasing as shown in Corollary 3.2, then
sup γκ
Z
t0 +amax
t0
v(αmax , t)p(ˆ u, t)dt ≤ γκamax
sup t0 ≤t≤t0 +amax
v(αmax , t)M (61)
Let supt0 ≤t≤t0 +amax v(αmax , t) = Vmax , then
sup γκ
Z
t0 +amax
t0
v(αmax , t)p(ˆ u, t)dt ≤ γκamax Vmax M. 17
(62)
ACCEPTED MANUSCRIPT
Divide (60) by γκamax M: (1 + C2 amax ) inf uˆ(0) − N γκamax M
Therefore, if Vmax ≤
(63)
CR IP T
Vmax ≤
(1 + C2 amax ) inf uˆ(0) − N , γκamax M 1 |C2 | < amax
then
0 < N < inf (u(a, t)) < inf (u(0, t)) .
t
AN US
a,t
Corollary 3.6. There exist an unique solution of u(a, t).
Proof. As stated in Corollary 3.3 the solution of u(a, t) is a system of ODEs with the form:
M
dˆ u = C2 uˆ − γκv(αmax , t)p(ˆ u, t) = F (ˆ u, t). dt
(64)
Let uˆ1 and uˆ2 be different, then:
CE
PT
ED
|F (ˆ u1 , t) − F (ˆ u2 , t)| = |C2 (ˆ u1 − uˆ2 ) − γκv(αmax , t)(p(ˆ u1 , t) − p(ˆ u2 , t))| . (65) R amax u ˆn Because p(ˆ un , t) = En (t) , where En (t) = a0 un (a, t)da. (65) further simplifies to: uˆ uˆ2 1 − . (66) |F (ˆ u1 , t) − F (ˆ u2 , t)| = C2 (ˆ u1 − uˆ2 ) − γκv(αmax , t) E1 (t) E2 (t)
AC
As a result, it follows that: u ˆ E (t) − u ˆ E (t) 1 2 2 1 |F (ˆ u1 , t) − F (ˆ u2 , t)| = C2 (ˆ u1 − uˆ2 ) − γκv(αmax , t) E1 (t)E2 (t) uˆ1 E2 (t) − uˆ2 E1 (t) ≤ |C2 ||ˆ u1 − uˆ2 | + γκv(αmax , t) E1 (t)E2 (t) uˆ E (t) − uˆ E (t) 1 2 2 1 ≤ |C2 ||ˆ u1 − uˆ2 | + |γκv(αmax , t)| E1 (t)E2 (t) (67) 18
ACCEPTED MANUSCRIPT
CR IP T
In the previous Corollary, we have proved that under suitable initial conditions and parameters, 0 < N < inf t u(a, t) < u(a, t) < supt u(a, t). Which implies that En (t) is bounded. Let A = inf t En (t) and B = supt En (t). Without loss of generality. assume E1 (t) ≤ E2 (t). It follows that (67) is bounded above by: supt E2 (t) |ˆ u1 − uˆ2 | inf t E1 (t) inf t E2 (t) t B ≤ |C2 ||ˆ u1 − uˆ2 | + sup |γκv(αmax , t)| 2 |ˆ u1 − uˆ2 | A B u1 − uˆ2 |. = |C2 | + sup |γκv(αmax , t)| 2 |ˆ A
AN US
sup |C2 ||ˆ u1 − uˆ2 | + sup |γκv(αmax , t)|
If we define:
B Ω = |C2 | + sup |γκv(αmax , t)| 2 , A then we have established that:
(69)
(70)
M
|F (ˆ u1 , t) − F (ˆ u2 , t)| ≤ Ω|ˆ u1 − uˆ2 |.
(68)
ED
We have demonstrated that each of the infinite number of ODEs composing the solution of u(a, t) is Lipschitz continuous. Applying the Picard-Lindel¨of theorem, we arrive at the conclusion that there exist an unique solution for each of the ODEs composing the solution of u(a, t) and therefore there exist an unique solution of u(a, t).
PT
Corollary 3.7. The unique solution for v(α, t) is stable for a set of parameters and initial conditions.
AC
CE
Proof. Assume that there are two solutions of the system v1 (α, t) and v2 (α, t) with different boundary conditions:v1 (0, t),v1 (α, 0),v2 (0, t) and v2 (α, 0). Recall Corollary 3.1 which states that the solution of v(α, t) is a system of ODEs with the following form: ˆ (t) = V t + a0 α (71) vˆ(t) = v(tˆ0 ) exp(C1 t) vˆ(t0 ) = v(a0 , 0). Using (41), it is clear that, v1 (α, t) = v1 (t − V t, 0) exp(C1 t) on its characteristic curves and v2 (α, t) = v2 (t − V t, 0) exp(C1 t) on its characteristic curves. 19
ACCEPTED MANUSCRIPT
Therefore, the difference between v1 (α, t) and v2 (α, t) at any given α and t can be expressed as the following if t < αmax : (72)
CR IP T
|v1 (α, t) − v2 (α, t)| = |(v1 (t − V t, 0) − v2 (t − V t, 0)) exp(C1 t)| .
If t > αmax , using (42), the difference between v1 (α, t) and v2 (α, t) at any given α and t can be expressed as the following: |v1 (α, t) − v2 (α, t)| = |(v1 (0, t0 ) − v2 (0, t0 )) exp(C1 (t − t0))| , where
(73)
α . (74) V Recall Corollary 3.2, which states that there exist a set of parameters such that v(α, t) is bounded. Therefore, according to (72) and (73) it is clear that the solution of v(α, t) depends continuously on its initial and boundary conditions, thus v(α, t) is stable.
AN US
t0 = t −
Corollary 3.8. The unique solution for u(a, t) is stable.
AC
CE
PT
ED
M
Proof. Assume that there are two solutions of the system u1 (a, t) and u2 (a, t) with different boundary conditions: u1 (0, t), u1 (a, 0), u2 (0, t) and u2 (a, 0). Then the difference between the total numbers of RBCs in u1 (a, t) and u2 (a, t) at an arbitrary time point t can be expressed in the following form as demon-
20
ACCEPTED MANUSCRIPT
amax
u1 (a, t1 ) − u2 (a, t1 ) da
a0
a0
+
Z
t2
t Z 1t2
CR IP T
strated in (1): Z Z amax u1 (a, t2 ) − u2 (a, t2 ) da =
u1 (a0 , t) − u2 (a0 , t) dt
u1 (amax , t) − u2 (amax , t) dt Z amax Z t2 v(αmax , t)r(a)p(u1 (a, t)) dt da −γ t1 a0 Z amax Z t2 +γ v(αmax , t)r(a)p(u2 (a, t)) dt da −
a0 t1 Z t2 Z T a max X i=1 Q
t1
amax Z t2
a0
i=1 amax Z t2
Z
a0
t1
t1
wi (t)θi u1 (a, t) − u2 (a, t) dt da
si (t)ψi u(a, t) − u2 (a, t) dt da
h(a)u1 (a, t) − u2 (a, t) dt da .
ED
−
a0
XZ
M
−
t1
AN US
−
AC
CE
PT
Let W = u1 (a, t) − u2 (a, t), then (75) simplifies to: Z amax Z amax Z t2 W (a, t2 ) da = W (a, t1 ) da + W (a0 , t) dt a0 a0 t1 Z t2 − W (amax , t) dt t1 Z amax Z t2 −γ v(αmax , t)r(a)p(u1 (a, t)) dt da a0 t1 Z amax Z t2 +γ v(αmax , t)r(a)p(u2 (a, t)) dt da −
a0 t1 Z Z t2 T a max X i=1
a0
t1
21
wi (t)θi W (a, t) dt da
(75)
ACCEPTED MANUSCRIPT
−
amax Z t2
si (t)ψi W (a, t) dt da
t1 i=1 a0 Z amax Z t2
h(a)W (a, t) dt da .
t1
a0
Assuming r(a) is a constant,then the following is true: Z t2 Z amax Z t2 v(αmax , t)r(a)dt. v(αmax , t)r(a)p(u1 (a, t)) dt da = γ γ a0
t1
t1
(76)
CR IP T
−
Q Z X
(77)
t1 T X Z amax Z t2 i=1 Q
a0
XZ
M
−
AN US
Substituting (77) into (76) we arrive at the following system: Z amax Z amax Z t2 W (a, t2 ) da = W (a, t1 ) da + W (a0 , t) dt a0 a0 t1 Z t2 − W (amax , t) dt
ED
−
−
wi (t)θi W (a, t) dt da
t1
amax Z t2
si (t)ψi W (a, t) dt da
t1 i=1 a0 Z amax Z t2
h(a)W (a, t) dt da .
a0
(78)
t1
CE
PT
If W (a, t) and its partial derivatives are continuous functions then (78) can be rearranged into an transport system as demonstrated in (36): ∂W ∂W + = C2 W (a, t), ∂t ∂a
(79)
AC
subject to the following initial and boundary condition: W (a, 0) =u1 (a, 0) − u2 (a, 0), W (0, t) =u1 (0, t) − u2 (0, t).
(80)
As shown in Corollary 3.1,the solution to (79) is a system of ODEs with the following form: ( da = 1, dt (81) ˆ dW ˆ, = C2 W dt
22
ACCEPTED MANUSCRIPT
which states that alone the characteristic curve given by The solution to (81) along the characteristic curve is:
da dW , dt dt
= exp(C2 t).
ˆ (t) = W ˆ (0) exp(C2 t), W
(82)
CR IP T
ˆ (0) is the initial condition for the ODE on the characteristic curve where W and lies either on W (a, 0) or W (0, a). Recall Corollary 3.5 and Corollary 3.4, which state that there exists a set of parameters such that u(a, t) is bounded. Therefore, according to (81) and (82), the solution to u(a, t) depends continuously on its initial and boundary conditions, thus u(a, t) is stable.
AN US
Lemma 3.9. The PDE system (26) is well-posed for a set of parameters and initial conditions.
M
Proof. Corollary 3.1 and Corollary 3.3 demonstrate that the solution of the two transport equation within (26) are systems of first order ODEs. Corollary 3.6 and Corollary 3.1 proved the existence and uniqueness of the solution to the first order ODEs composing the solution to the two transport equation within our system. Additionally, Corollary 3.7 and Corollary 3.8 demonstrate that the solution to (26) depends continuously on the initial and boundary conditions. Therefore, the PDE system (26) is well-posed.
ED
4. Numerical Solution
CE
where
PT
Recall that our system is a transport system with the form:
m(a, t) =
T X i=1
∂u ∂u + = −m(a, t)u(a, t), ∂t ∂a
wi (t)θi +
Q X i=1
γv(αmax , t)r(a) si (t)ψi + h(a) + R amax u(a, t) a0
(83) !
.
(84)
AC
Using the first order forward finite difference scheme at u(a, t), where the time interval is divided into N steps of size k (indicated with the superindex j), and the age interval is divided into M segments of size h (indicated with the subindex i). it is clear that: ∂u uj+1 − uji = i , ∂t k 23
(85)
ACCEPTED MANUSCRIPT
and
AN US
CR IP T
uj − uji−1 ∂u = i . (86) ∂a h Substituting the above Equation into (83), and let r = k/h, the system becomes: uji − uji−1 uj+1 − uji i + = −mji uji , (87) k h resulting in the following difference equation: j+1 j j i j ui = −kmj ui + (1 − r)ui + rui−1 , (88) u0i = g(i), j u0 = f (j). Using the above schema one can iteratively approximate all u(a, t) at t = k then t = 2k etc.
Lemma 4.1. j The above finite difference schema is numerically stable if r < 1 and max{ mi : ∀i, j}k < (2 − 2r).
PT
ED
M
Proof. Let uˆji be the numerical solution of (88). Then the error eji = uˆji − uji satisfies: j+1 j j i j ei = −kmj ei + (1 − r)ei + rei−1 , (89) e0i = 0, j e0 = 0. Let ej = max{ eji : ∀i} for each j ≥ 0. Then
AC
CE
j+1 e = (1 − r + kmij ) + rej i i−1 (90) i ≤ ( 1 − r + kmj + r)ej Thus ej+1 ≤ ( 1 − r − kmij + r)ej for all j ≥ 0. Let m ˆ = max{ mji : ∀i, j} for each i ≥ 0 and j ≥ 0. For r < 1 and k m ˆ < (2 − 2r), we have ej+1 ≤ ej for all j ≥ 0. This means that the error ej is not increasing when j is increasing. Therefore, (88) is stable. The same schema is also applied to v(α, t). Recall that: ∂v ∂v +V = −n(α, t)v(α, t), ∂t ∂α 24
(91)
ACCEPTED MANUSCRIPT
where n(α, t) =
T X
wi (t)φi +
Q X i=1
i=1
!
si (t)bi (t) .
(92)
CR IP T
Using the same finite difference schema one can iteratively approximate all v(α, t) at t = k then t = 2k etc. j+1 j j i j vi = −knj vi + (1 − rV )vi + rV vi−1 (93) vi0 = c(i) j v0 = d(j)
AN US
Which is also stable for rV < 1 and max{nij }k < (2 − 2rV ).
M
Additionally, using Taylor expansion, the value of w(t) and s(t) are approximated as following: dwi wi (t + k) = w(t) + k dt (t), R αmax dwi (94) v(a, t)dα, (t) = o (t) − β w (t) − τ w i i i i i dt α0 wi (0) = A,
PT
5. Results
dsi si (t + k) = s(t) + k dt (t), Rα dsi (t) = li (t) − δi si (t) − ϑi si α0max v(a, t)dα, dt si (0) = B.
(95)
ED
and
AC
CE
5.1. RBC Equilibrium Condition The erythropoietic process responds to events altering the concentration of RBCs and ensures that the number and age distribution of RBCs return to normal level upon perturbation. Figure 1 shows that using an arbitrary initial age distribution of RBCs, the system eventually returns to the basal RBC age distribution overtime. This result demonstrates that our system reflects the actual biological property of erythropoietic control, where the host responds to a perturbation in concentration of RBCs by altering accordingly, the rate of RBC production in the bone marrow, and the rate of RBCs entering circulation. This regulatory effect can be seen with the fluctuation of u(0, t) in Figure 1, which first responds to the lower than normal number of RBCs with a increase of RBCs entering circulation after a delay of 20 days. 25
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Figure 1: Using an arbitrary initial condition, the system returns to Equilibrium condition after 80 days.
AC
CE
PT
ED
5.2. Susceptibility As stated in Corollary 3.2, if the immune strength is weak, eventually the iRBC population grows without bounds and RBCs are depleted, resulting in severe anemia (Figure 2). Depletion of RBCs is characteristic of first time malaria infection which without treatment could result in mortality of the host. Within our model, this behavior can be achieved under at least two parameter sets. The first one is where C1 , the combined effector strength upon saturation of both innate and adaptive immune cells is weak, which allows for the depletion of RBCs. The second set is where C1 is large enough to halt the exponential growth of iRBCs, but the parameters controlling the speed at which innate immune cells, adaptive immune cells and adaptive immune effector strength growth are small, resulting in the depletion of RBCs before C1 can reach its maximum value. The second scenario reflects the case where a host is susceptible to severe anemia during primary infection, but during repeated infections, C1 had enough time to reach maximum value, and the host is resilient to severe anemia [3].
26
,
CR IP T
ACCEPTED MANUSCRIPT
AN US
Figure 2: The iRBC population grows exponentially due to weak immune strength
ED
M
5.3. Resilience
ln 1/γκ αmax /V
PT
Figure 3: iRBC and RBC co-exist when total immune strength is close to
AC
CE
1/γκ , the iRBC and RBC If the immune effector strength is close to αlnmax /V population coexist for a long period of time as shown in Figure 3. This state of resilience is characteristic of repeated malaria infection, where the host has already built up an immune strength against the parasite. However, sustained periodic coexistence is almost never observed in the real world, which is expected due to the fixed maximum immune strength required for sustained coexistence.
5.4. Resistance As stated in Corollary 3.2, if the immune strength reaches a greater value 1/γκ than αlnmax , then iRBC population decays constantly until the complete /V 27
CR IP T
ACCEPTED MANUSCRIPT
ln 1/γκ αmax /V
AN US
Figure 4: iRBC and RBC co-exist when total immune strength is greater than
removal of iRBC, as seen in Figure 4. This state of resistance is rarely observed without medical intervention during first infection; however, it is feasible to develop an immune response against malaria upon continuous and repeated exposure high enough to clear the parasite [16].
M
6. Discussion
AC
CE
PT
ED
In this paper we present a theoretical model to capture the complex interaction between iRBCs, RBCs and the host immune system during the course of a Plasmodium infection. Despite not knowing many of the parameters regarding the immune system’s effect on RBCs and iRBCs, our model was able to capture the general behavior of two major outcomes of infection with Plasmodium: (i) the depletion of host RBCs in a naive individual and (ii) the coexistence of RBCs, iRBCs, and the immune system in individuals with previous exposure. Our model states that the immune effector has to be equal to an exact number for sustained coexistence of the system, which on first look seems rather impossible as biological systems are rarely precise. This further coincides with known malaria infection dynamics where the oscillation of iRBC populations is rarely sustained or regular. To produce sustained oscillation over weeks or months, the immune effector strength only needs to lie in a range rather than being fixed. It is worth mentioning, that one of the major simplifications of the model was the assumption of r(a), the infection probability of RBCs of different age, to be constant. For sufficiently regular r(a), the well-posedness results can be derived in a similar fashion. However, the rigid trichotomy which 28
ACCEPTED MANUSCRIPT
AC
CE
PT
ED
M
AN US
CR IP T
conveniently characterizes the system is not inherited. The simplified system in (28) results in uncoupled dynamics for iRBCs. This could imply that infected cells may appear even in the absence of uninfected cells; however, this condition would not be physiological. Given realistic initial and boundary conditions, the simplified system permits analysis and approximates well the reality of a malaria infection. We plan to calibrate our model in the future using time course data from the Malaria Host-Pathogen Interaction Center (MaHPIC), which is a large collaborative project between Emory University, The University of Georgia, The Georgia Institute of Technology, and the Centers for Disease Control and Prevention. The MaHPIC project aims to obtain time courses of -omic data (transcriptomic, metabolomic, lipidomic, proteomic) from non-human primates infected with several different species of Plasmodium parasites during multiple experiments, with initial experiments lasting for about 100 days. Using data produced by the MaHPIC group, which include daily RBC counts, iRBC counts and immune cell counts over the 100-day period, we plan to calibrate the full model, including the different infection rates of different age groups of RBCs and the different effector strengths against different age groups of iRBCs. We expect to find discrepancies between our model prediction and the actual experimental data; such discrepancy will not only provide biological insight regarding the hematodynamics of malaria infection, but will also allow us to re-iterate our model. Our current PDE model is by no means a complete reflection of the complex within-host dynamics of malaria due to its vast simplification of the immune system’s effect on the iRBC and RBC populations along with its complete disregard of the within immune system interaction. However, the model could be expanded to also contain these features and serve as a liaison where the molecular features of the system can be reflected at the cellular level. We plan to use the molecular level data generated by MaHPIC to further expand our model to incorporate the more subtle aspects of the within-host interaction of malaria infection. Acknowledgments The authors would like to thank Jacob A. Grey for his help in refining the corollaries 3.4 through 3.7, and also Stacey Lapp, Tracey Lamb, Rabindra Tirouvanziam, Chet Joyner, and Esmeralda Meyer for fruitful discussions. Two anonymous reviewers helped to improve this manuscript significantly. 29
ACCEPTED MANUSCRIPT
CR IP T
This research has been supported Malaria Host-Pathogen Interaction Center - MaHPIC, NIH’s NIAID contract HHSN272201200031C, and the International Centers for Excellence in Malaria Research - Center for non-Amazonian regions of Latin America (CLAIM), NIH’s NIAID cooperative agreement U19AI089702. References
[1] Z. Bozdech, M. Llin´as, B. L. Pulliam, E. D. Wong, J. Zhu, J. L. DeRisi, The transcriptome of the intraerythrocytic developmental cycle of plasmodium falciparum, PLoS biology 1 (1) (2003) e5.
AN US
[2] W. E. Collins, M. Warren, J. S. Sullivan, G. G. Galland, Plasmodium coatneyi: observations on periodicity, mosquito infection, and transmission to macaca mulatta monkeys., The American journal of tropical medicine and hygiene 64 (3) (2001) 101–110.
M
[3] A. Moreno, M. Cabrera-Mora, A. Garcia, J. Orkin, E. Strobert, J. W. Barnwell, M. R. Galinski, Plasmodium coatneyi in rhesus macaques replicates the multisystemic dysfunction of severe malaria in humans, Infection and immunity 81 (6) (2013) 1889–1904.
ED
[4] I. M. Rouzine, F. E. McKenzie, Link between immune response and parasite synchronization in malaria, Proceedings of the National Academy of Sciences 100 (6) (2003) 3473–3478.
PT
[5] D. A. Moreno-P´erez, J. A. Ru´ız, M. A. Patarroyo, Reticulocytes: Plasmodium vivax target cells, Biology of the Cell 105 (6) (2013) 251–260.
CE
[6] D. J. Perkins, T. Were, G. C. Davenport, P. Kempaiah, J. B. Hittner, J. M. Ong’echa, Severe malarial anemia: innate immunity and pathogenesis, International journal of biological sciences 7 (9) (2011) 1427.
AC
[7] T. J. Lamb, J. Langhorne, The severity of malarial anaemia in plasmodium chabaudi infections of balb/c mice is determined independently of the number of circulating parasites, Malaria journal 7 (1) (2008) 68. [8] N. M. Douglas, N. M. Anstey, P. A. Buffet, J. R. Poespoprodjo, T. W. Yeo, N. J. White, R. N. Price, et al., The anaemia of plasmodium vivax malaria, Malar J 11 (135.10) (2012) 1186. 30
ACCEPTED MANUSCRIPT
[9] M. Mackey, J. Milton, Feedback, delays and the origin of blood cell dynamics, University of Minnesota. Institute for Mathematics and Its Applications, 1990.
CR IP T
[10] J. Belair, M. C. Mackey, J. M. Mahaffy, Age-structured and two-delay models for erythropoiesis, Mathematical Biosciences 128 (1) (1995) 317– 346. [11] J. M. Mahaffy, J. Belair, M. C. Mackey, Hematopoietic model with moving boundary condition and state dependent delay] applications in erythropoiesis, Journal of Theoretical Biology 190 (2) (1998) 135–146.
AN US
[12] N. J. Savill, W. Chadwick, S. E. Reece, Quantitative analysis of mechanisms that govern red blood cell age structure and dynamics during anaemia, PLoS computational biology 5 (6) (2009) e1000416. [13] J. J. Thibodeaux, Modeling erythropoiesis subject to malaria infection, Mathematical biosciences 225 (1) (2010) 59–67.
M
[14] A. S. Ackleh, B. Ma, J. J. Thibodeaux, A second-order high resolution finite difference scheme for a structured erythropoiesis model subject to malaria infection, Mathematical biosciences.
PT
ED
[15] C. Metcalf, A. Graham, S. Huijben, V. Barclay, G. Long, B. Grenfell, A. Read, O. Bjørnstad, Partitioning regulatory mechanisms of withinhost malaria dynamics using the effective propagation number, Science 333 (6045) (2011) 984–988.
AC
CE
[16] D. L. Doolan, C. Dobano, J. K. Baird, Acquired immunity to malaria, Clinical microbiology reviews 22 (1) (2009) 13–36.
31
ACCEPTED MANUSCRIPT
APENDIX: Table of Parameters and Variables Description
si (t)
h(a) f (t, ϕt ) r(a)
AC
CE
PT
ED
li (t)
M
p(u(a, t))
oi (t)
Source
n/a n/a n/a
n/a n/a n/a
n/a n/a n/a n/a
n/a n/a n/a n/a
n/a
n/a
0.02 0.03
MaHPIC
AN US
u(a, t) v(α, t) ϕ(t) wi (t)
bi (t)
Value
CR IP T
t a α
Unit Independent Variables Time. Day Age of RBC. Day Age of iRBC. 10− 1Day Dependent Variables Circulating RBC age distribution. cells/ul/day iRBC age distribution. cells/ul/day Circulating RBC concentration. cells/ul Circulating innate immune cell concencells/ul tration. Circulating adaptive immune cell con- cells/ul centration. Functions Percentage of RBC that leaves circula1/day tion per day, usually a constant. Rate at which RBC enters circulation. cells/ul/day Success rate of merozoite invading 1/day RBC of age. a Probability of a merozoite infecting a dimensionless RBC of age. a at time t Rate at which adaptive immunity cell 1/(day · cells/ul) i destroys iRBC. Rate at which innate immune cell encells/ul/day ters circulation. Rate at which adaptive immune cell cells/ul/day enters circulation.
32
n/a [0, 1]
n/a Literature
n/a
n/a
n/a
Estimated
n/a
Estimated
n/a
Estimated
ACCEPTED MANUSCRIPT
ψi φi βi δi τi ϑi λi $i σi i %i νi ηi ωi Mi
CE
Ri
xii
AC
ε ς
Td
Value
Source
10
n/a
6400000
MaHPIC
25
Literature
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
CR IP T
θi
M
γ
ED
ϕ0
PT
V
Unit Parameters Speed at which parasite ages, equal to 1/(day) 1/α. Normal circulating RBC concentracells/ul tion, equal to ϕ(0). Average number of merozoite pro- dimensionless duced by a single iRBC. Rate at which innate immunity cell i 1/(day · cells/ul) destroys RBC. Rate at which adaptive immunity cell 1/(day · cells/ul) i destroys RBC. Rate at which innate immunity cell i 1/(day · cells/ul) destroys iRBC. Rate at which innate immunity cell de- 1/day cays. Rate at which adaptive immunity cell 1/day decays. Rate at which innate immunity cell de- 1/(day · cells/ul) cay due to contact with iRBC. Rate at which adaptive immunity cell 1/(day · cells/ul) decay due to contact with iRBC. Coefficient for change of adaptive im- 1/cells/ul munity effector strength. Normal rate of production of innate cells/ul/day immune cell i. Normal rate of production of adaptive cells/ul/day immune cell i. Maximum rate of production of innate cells/ul/day immune cell i. Maximum rate of production of adapcells/ul/day tive immune cell i. Maximum adaptive immunity effector 1/day/cells/ul i strength. Coefficient for change of production of 1/cells/ul innate immune cell i. Coefficient for change of production of 1/cells/ul adaptive immune cell i. Amount of parasite where the produc- cells/ul tion of innate immune cell i. increase the most Amount of parasite where the produc- cells/ul tion of adaptive immune cell i. increase the most Amount of parasite where the produc- cells/ul tion of adaptive immune cell i. effector strength increase the most Coefficient for change of erythro- 1/cells/ul poiesis. Number of RBC entering blood cells/ul/day stream. Delay of hematopoietic response. day
AN US
Description
33
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
n/a
Estimated
8000
MaHPIC
10 20
Literature
