Interactions Between the Frequency of the Duffy Antigen and the Dynamics of P. vivax Malaria Infections

Interactions Between the Frequency of the Duffy Antigen and the Dynamics of P. vivax Malaria Infections Elizabeth Ghartey1,a,*, Gautam Rai2,a,*, Dasha Selivonenko3,a,*, Rachel K. Wissenbach (Mentor)2,a, Lucero Rodriguez Rodriguez (Mentor)2,a, and Joan Ponce (Mentor)2,a *These authors contributed equally to the work. 1University of Arizona, Tucson AZ, United States 2Arizona State University, Tempe AZ, United States 3Rice University, Houston TX, United States aSimon A. Levin Mathematical, Computational, and Modeling Sciences Center: Quantitative Research for the Life & Social Sciences Program, Arizona State University, Tempe AZ, United States July 19, 2024 Abstract The malarial parasite Plasmodium vivax has infected humans for millennia. As such, alleles like the FY*O allele in the Duffy Antigen Receptor for Chemokines (DARC) and the sickle cell allele (HbS) have been naturally selected for in malaria-endemic regions because they confer resistance to malaria, thereby increasing the survival and reproductive success of individuals carrying these protective alleles. As resistance becomes more common in a population, malaria incidence is expected to decline, reducing the evolutionary pressure for additional resistance. In this work, we explore the interaction between these processes. We construct a model with seasonality that tracks the frequency of Duffy genotypes and leverages fast/slow dynamics to analyze the coupled dynamics of malaria transmission and changes in the gene frequency of the DARC genotype. Specifically, we investigate how the burden of malaria changes with the fractions of people with the various DARC genotypes. We derive the basic reproduction number as the threshold condition for the stability of the disease-free equilibrium and interpret the R0as a weighted sum of the cases generated by infected individuals of each genotype. Additionally, we calibrate our model using data from the Amazonas region in Brazil, which has a polymorphic population with respect to DARC, and still reports a substantial number of P. vivax cases. Analysis of our model determines the proportion of the population that must be Duffy-negative in order for the entire population to be protected against P. vivax without any further interventions. Furthermore, we assess how different proportions of Duffy-negative individuals influence the monthly incidence of P. vivax cases. 1 Contents 1 Introduction 3 1.1 Influence of Genetics on the Epidemiology of Malaria ............. 3 1.2 Previous Mathematical Models of Natural Protection Against Malaria .... 4 2 Model 4 2.1 Model Formulation ................................ 4 2.2 Human population genetic model ........................ 8 2.3 Re-scaled model .................................. 10 2.4 Slow Dynamical System ............................. 10 2.4.1 Manifold of Slow System ......................... 11 2.4.2 Determining Fitness of the X-gene ................... 11 2.5 Fast Dynamical System .............................. 13 2.6 The Next Generation Matrix and the Basic Reproduction Number of the Fast System ....................................... 14 2.7 Conditions for Stability of Disease-Free Equilibrium .............. 15 2.8 Location of Nontrivial Equilibrium ....................... 16 2.9 Conditions for Stability of Endemic Equilibrium ................ 18 3 Sensitivity Analysis 20 3.1 Calculating The Elasticity of R0......................... 21 3.2 Global Sensitivity on Fast System with Biting Rate Seasonality ....... 21 4 Parameter Estimation 22 4.1 Relative Risk of Infection Based on Genotype ................. 22 4.2 Fitting the Fast System with Seasonality .................... 24 5 Results 26 6 Discussion 29 7 Acknowledgements 30 8 Appendix 31 8.1 Mosquito Model Re-scale (3) ........................... 31 2 1 Introduction 1.1 Influence of Genetics on the Epidemiology of Malaria Malaria is an infectious disease spread by mosquitoes that are infected with Plasmodium parasites. Whereas P. falciparum is the most dangerous strain and the most prevalent in Africa, P. vivax is the most geographically widespread strain of malaria and the most common in the rest of the world [13]. Malaria is a perennial affliction to many tropical countries: in 2022, it caused an estimated 249 million cases globally, resulting in more than six hundred thousand deaths and billions of dollars in global costs [23]. The selective pressure induced by malaria has demonstrably influenced the genetic composition of tropical populations. Examples include hemoglobinopathies like the sickle cell trait in Sub-Saharan Africa and β-thalassemia in the Eastern Mediterranean. The sickle cell trait is caused by a missense mutation of the β-globin gene, resulting in the abnormal synthesis of the β-globin protein. The vast majority of people with sickle cell trait and sickle cell disease live in Sub-Saharan Africa. Similarly, β-thalassemia is characterized by obstructed synthesis of the β-globin chain, which causes numerous phenotypes that may be asymptomatic or debilitating. People with sickle cell trait and some phenotypes of β-thalassemia have significant protection against P. falciparum malaria compared to people without either condition [6,8, 11]. Another example of an adaptation against malaria is the absence of the Duffy antigen receptor for chemokines (DARC) on human red blood cells. P. vivax protozoa rely on the interaction between DARC and the Plasmodium vivax Duffy-binding protein (PvDBP) to invade human reticulocytes. Having the Duffy-negative trait thereby removes a mode of invasion into reticulocytes by P. vivax protozoa. Nonetheless, infections by P. vivax have been observed in Duffy-negative individuals, but the mechanism of infection in these cases are not fully understood [20]. The transmission patterns of malaria are highly influenced by the Duffy-negative trait. Duffy-negativity is very common in sub-Saharan Africa, so the fraction of the population at risk of vivax malaria is very small. Consequently, the majority of malaria infections are caused by P. falciparum. Although P. vivax infections is exceedingly rare in sub-Saharan Africa, it is the most common malarial parasite in the rest of the world [10]. Our system models the dynamics of P. vivax malaria in the Amazonas state of Brazil, because this population is much more polymorphic with respect to the Duffy antigen trait than most parts of the world. Kano et al. [16] conducted a longitudinal study on the influence of DARC on people’s susceptibility to vivax malaria in Rio Pardo, which is located in the Brazilian state of Amazonas. In this study, they sampled antibody response to vivax malaria, and disaggregated it by Duffy genotype to ascertain the influence of natural immunity. We used data from the longitudinal work of Kano et al to determine the populations of interest in our study. 3 1.2 Previous Mathematical Models of Natural Protection Against Malaria Previous mathematical models of malaria transmission dynamics frequently use epidemiological compartmental models like the SIR model. In this paper, we use ordinary differential equations to characterize the dynamics of vivax malaria in Rio Pardo, Brazil, with relation to the Duffy antigen. We reference previous population genetics models that relate to the sickle cell trait. Feng et al. [7] used a system of ordinary differential equations to describe the relationship between malaria transmission dynamics and the genetic frequency of the S-gene that causes sickle cell anemia. Furthermore, Feng and Castillo-Chavez [6] used an expanded system of ordinary differential equations with three elements to characterize the complex dynamics of the S-gene beyond the standard two-element system. We use the Rio Pardo area studied by Kano et al. [16], because it is more polymorphic than most parts of the world. This has implications for transmission dynamics of vivax malaria in the region. The paper is structured as follows: we introduce the terms, parameters, and variables of our mathematical model. We then describe our model and its dynamics. We formulate the slow system, from which we deduct its manifold and the fitness of the X-gene. We use this information to simulate the fast system and to determine its next generation matrix, basic reproduction number, and equilibria. We conduct a global sensitivity analysis on the model, eventually including biting rate seasonality, and then estimate our parameters. We then conclude the paper by assessing the results and implications of our model. 2 Model 2.1 Model Formulation There are four possible phenotypes of the Duffy antigen: Duffy-positive, where two of the alleles (FY*A and/or FY*B) alleles are present; the two heterozygous Duffy-positive traits, where one of the two alleles is present and the other is silent; and the homozygous Duffynegative trait, where neither the FY*A or FY*B allele is expressed. For our paper, we use the following nomenclature for the Duffy blood group [12,16,18]: Meaning Phenotype Alleles & Genotype Classification Duffy-Positive Fy(a+b+) FY*A/FY*B XX Fy(a+b−) FY*A/FY*A XX FY*A/FY*O XO Fy(a−b+) FY*B/FY*B XX FY*B/FY*O XO Duffy-Negative Fy(a−b−) FY*O/FY*O OO In this paper, we use the term XX to refer to the presence of two expressed Duffy alleles, XO for the presence of one expressed Duffy allele and one silent Duffy allele, and OO for 4 the presence of two silent Duffy alleles. The subscript idenotes the genotype: i= 1 for XX, i= 2 for XO, and i= 3 for OO. We structure our model similarly to the one that Feng and Castillo-Chavez [6] used to describe the dynamics of the sickle cell trait. Distinct differences between our models exist: whereas the sickle cell model developed by Feng incorporates a death rate associated with the sickle cell trait, our model does not implement a death rate due to the Duffy antigen. Furthermore, our model considers the impact of seasonality on our numerical simulations. The variables of our initial model are described in Table 1 and a summary of the parameters is presented in Table 2.1. Table 1: Variables and description for original model. Our values for Piand wiare based in Rio Pardo, Brazil. State Variable Description uiUninfected individuals of genotype i viInfected individuals of genotype i zFraction of infected mosquitoes xiFraction of uninfected individuals of genotype i yiFraction of infected individuals of genotype i wiFraction of individuals in each genotype NTotal population SmNumber of uninfected mosquitoes ImNumber of infected mosquitoes MTotal number of mosquitoes 5 Table 2: Parameters, descriptions, values, units, and source. The values for Piand wiare for Rio Pardo, Brazil. Human Parameters Description Value Unit Source P1 Fraction of total births with XX genotype 0.6657 people [16] P2 Fraction of total births with XO genotype 0.3003 people [16] P3 Fraction of total births with OO genotype 0.0339 people [16] w1Fraction of XX people 0.6623 N/A [16] w2Fraction of XO people 0.3072 N/A [16] w3Fraction of OO people 0.0304 N/A [16] NTotal Population 690 People [16] bBase human birth rate 1.26e-5 1/day [4] KDensity dependent reduction in birth rate 10000 people [6] θ1 Probability of infection from infected bite for XX 0.78 N/A [5] θ2 Probability of infection from infected bite for XO 0.78 N/A [5] θ3 Probability of infection from infected bite for OO 0.00702 N/A [5,16] γ1Recovery Rate for XX 1/24 1/day [27] γ2Recovery Rate for XO 1/24 1/day [27] γ3Recovery Rate for OO 1/24 1/day [27] mNatural human death rate 3.78e-5 1/day [4] α1 Disease-induced death rate for XX People 0.0003 Deaths/infection [19] α2 Disease-induced death rate for XO People 0.0003 Deaths/infection [19] α3 Disease-induced death rate for OO People 0.0003 Deaths/infection [19] Mosquito Parameters aBiting Rate 0.5 1/day [24] cMosquitoes per Human 2.3 - 4.36 Mosquitoes/Human [25] ϕ1 Probability a mosquito is infected from XX 0.326 N/A [1] ϕ2 Probability a mosquito is infected from XO 0.326 N/A [1] ϕ3 Probability a mosquito is infected from OO 0.30 N/A [22] δMosquito Death Rate 0.315 1/day [2] 6 Let u(t) be the total number of uninfected people and v(t) be the total number of infected people at time t. Thus, the total population is given by N(t) = u(t) + v(t). Furthermore, let the subscript i= 1 describe individuals with the Duffy genotype XX; i= 2 describe individuals with the XO genotype; and i= 3 describe individuals with the OO genotype. For example, u1denotes the number of uninfected people whose Duffy genotype is XX, and v2denotes the number of infected people with the XO genotype. Lastly, let Sm(t) be the total number of susceptible mosquitoes and Im(t) be the total number of infected mosquitoes, where the total number of mosquitoes in the population is given by M(t)≡Sm(t) + Im(t) at time t. The mathematical model is derived based on the following assumptions: 1. New uninfected people of genotype iare born at the rate PiNb(N), where b(N) represents the per capita birth rate as a function of the total population, and Piis the fraction of births with genotype i. Arbitrarily, we define b(N) as b(1 −N K), where b is the base birth rate and Kis the carrying capacity. We assume that people die of non-disease-induced causes at the rate of m. Infection occurs in humans at the rate of aθic, where cis the mosquito-to-human ratio, ais the mosquito biting rate and θiis the probability a human becomes infected after being bitten by an infected mosquito. Lastly, infected individuals recover and enter the uninfected class at the rate γi. 2. Individuals enter the infected compartment at a rate aθic. We assume individuals either die of natural causes at a rate m, or die due to malaria at a rate αi. Additionally, infected individuals may recover and leave the infected class at a rate γi. 3. The mosquitoes are born into the susceptible class at rate bm. The term SmP3 i=1 haϕi vi Ni describes the transition from susceptible to infected mosquitoes, where ais the mosquito biting rate, and ϕirepresents the probability that a mosquito becomes infected after biting an infected human of genotype i. Note that vi Ncan be thought of as the probability that a mosquito bites a person of genotype igiven that they bite an individual. Lastly, we assume both susceptible and infected mosquitoes die at the rate of δ. We then add the susceptible and infected mosquito classes to create the compartment dM dt=bm−δM, which we use to monitor the total population of mosquitoes. 7 Thus we have the system: dui dt=PiNb(N) | {z } Birth −mui |{z} Natural death −aθiczui | {z } Infection +γivi |{z} Recovery dvi dt=aθiczui | {z } New infections −mvi |{z} Natural death −γivi |{z} Recovery −αivi |{z} Disease-induced death dSm dt=bm−Smaϕ1 v1 N(t)+aϕ2 v2 N(t)+aϕ3 v3 N(t) | {z } Rate of infection by blood meal −δz |{z} Natural death dIm dt=Smaϕ1 v1 N(t)+aϕ2 v2 N(t)+aϕ3 v3 N(t) | {z } Rate of infection by blood meal −δIm |{z} Natural death (1) The flowchart below describes how people and mosquitoes move between their respective classes. Birth Figure 1: Flowchart of the population dynamics of humans (above) and mosquitoes (below), as described in (3). 2.2 Human population genetic model Our model requires distinguishing between the dynamics of each population’s genotype, because malaria dynamics occur on two time scales. Thus, we used Punnett squares to determine the fraction of total births of each genotype, in the Rio Pardo locality of Amazonas, Brazil, denoted by Pi. We combine the expressed alleles FY*A and FY*B into a singular theoretical allele X. Likewise, we use O to represent the silent allele FYES (where ES means “erythrocyte silent”). Let wibe the fraction of the total population by genotype. An example calculation of P1is below. Note that since P1is the fraction of births with genotype XX, we only consider the possible crosses that can produce a child with XX when calculating P1. 8 Parent 2 (w1) Parent 1 (w1) X X X XX XX X XX XX =⇒ w1xw1 Genotype Probability XX 1 Parent 2 (w1) Parent 1 (w2) X X X XX XX O XO XO =⇒ w1xw2 Genotype Probability XX 1/2 XO 1/2 Parent 2 (w2) Parent 1 (w1) X O X XX XO X XX XO =⇒ w2xw1 Genotype Probability XX 1/2 XO 1/2 Parent 2 (w2) Parent 1 (w2) X O X XX XO O XO OO =⇒ w2xw2 Genotype Probability XX 1/4 XO 1/2 OO 1/4 P1= 1 (w1w1) + 1 2(w1w2) + 1 2(w2w1) + 1 4(w2w2) = (w1)2+ (w1w2) + 1 4(w2)2. The other formulations of Piare determined using a similar method. Thus, the formulas are as follows: P1= (w1)2+1 2(2w1w2) + 1 4(w2)2, P2=1 2(2w1w2) + 1 2(2w2w3) + 1 2(w2)2+ (2w1w3), P3= (w3)2+1 4(w2)2+1 2(2w2w3). To simplify our analysis, we let w3= 1 −w1−w2. Then, our Piexpressions are reduced: P1= (w1)2+ (w1w2) + 1 4(w2)2, P2= (w1w2)+(w2w3) + 2(w1w3) + 1 2(w2)2 P3= (w3)2+ (w2w3) + 1 4(w2)2 (2) The assumptions above lead to the following nonlinear ordinary differential equations, which model the disease dynamics of each genotype: 9 It can be shown that the matrix MD−1has two eigenvalues of 0 and the remaining two are λ1, λ2=±sw1 βh1βv1 δγ3 +w2 βh2βv2 δγ2 +w3 βh3βv3 δγ3 . Making the substitution R0=w1βh1βv1 δγ1+w2βh2βv2 δγ2+w3βh3βv3 δγ3yields λ1, λ2=±pR0. Note that λ2=−√R0is always less than 1 since R0is always greater than or equal to 0. Additionally, λ1<1 when R0<1 and λ1>1 when R0>1. By Theorem 2, the real part of all the eigenvalues of J(E0) are negative when R<1, meaning the disease-free equilibrium is locally asymptotically stable when R0<1. Likewise, E0is unstable when R0>1. 2.8 Location of Nontrivial Equilibrium The nontrivial equilibrium E∗= (y∗ 1, y∗ 2, y∗ 3, z∗) is a solution to the following system. 0 = βh1z∗(w1−y∗ 1)−γ1y∗ 1 0 = βh2z∗(w2−y∗ 2)−γ2y∗ 2 0 = βh3z∗(w3−y∗ 3)−γ3y∗ 3 0 = (1 −z∗)(βv1y∗ 1+βv2y∗ 2+βv3y∗ 3)−δz∗ (20) For convenience, we define Thi =βhi γi . Substituting and simplifying yields the following expression for y∗ i: y∗ i=βhiwiz∗ γi+βhiz∗=Thiwiz∗ Thiz∗+ 1.(21) To find the solutions for z∗, we substitute out y∗ 1, y∗ 2,and y∗ 3in (20) and divide out the trivial solution of z∗= 0. Additionally, we let w=w1+w2and assume Th1=Th2and R01 =R02. This yields the quadratic: 0 = k2(z∗)2+k1(z∗) + k0, where k2=Th1Th3+R03Th1(1 −w) + R01Th3w, k1=Th3+Th1+R03(1 −Th1)(1 −w) + R01w(1 −Th3), k0= 1 −R03(1 −w)−R01w. (22) Note that k0can be rewritten as k3= 1 −R0, and let h(z) = k2(z)2+k1(z) + k0. Note that z∗is by definition a solution to the equation h(z) = 0. Lemma 1. If 0<R0<1then k1>0. 16 Proof. Let 0 <R0<1 and define T0= max(Th1, Th3). Based on these two statements, Th3+Th1−T0R0≥0.(23) Additionally, since T0≥Th1and T0≥Th3, we know R03T0(1 −w) + R01T0w≥R03Th1(1 −w) + R01Th3w. (24) From our definition of k1(22), k1=Th3+Th1+R03(1 −Th1)(1 −w) + R01w(1 −Th3), =Th3+Th1+R03(1 −w) + R01w−R03Th1(1 −w)−R01Th3w, =Th3+Th1+R0−R03Th1(1 −w)−R01Th3w. From here, using relation (24) yields k1≥Th3+Th1+R0−R03T0(1 −w)−R01T0w, k1≥Th3+Th1+R0−T0R0, and relation (23) implies Th3+Th1−T0R0+R0≥R0. Combining the two results above yields k1≥R0. Therefore, if 0 <R0<1, then k1>0. Theorem 3. Assuming all parameters are biologically possible, the fixed point E0in the fast system (12a) is biologically impossible when R0<1and biologically possible when R0>1. Proof of Theorem 3.Let all parameters be biologically possible and 0 <R0<1. From Equation 22 it follows that k2>0 and k0>0, due to the respective conditions above. Additionally, Lemma 1 implies, k1>0. Let z∗ 1, z∗ 2be the two solutions to the quadratic h(z) = 0. Then, by Vieta’s formulas, the solutions to h(z) = 0 satisfy z∗ 1·z∗ 2=k0 k2 z∗ 1+z∗ 2=−k1 k2 and using the conditions k0, k1, k2>0, they simplify to z∗ 1·z∗ 2>0z∗ 1+z∗ 2<0 The inequality on the left rules out the possibility that exactly one of the roots of h(z) is positive and the inequality on the right implies that the roots cannot both be positive. Therefore, neither z∗ 1or z∗ 2are positive. As such, when 0 <R0<1, the nontrivial equilibrium is not biologically possible since E∗would require a negative or complex fraction of mosquitoes to be infected, neither of which are possible in the real world. To prove the second part of Theorem 3, let all parameters be biologically possible and let R0>1. Similar to earlier, these two conditions imply k2>0, k1<0 and k0<0.Similar to before, applying these conditions to Vieta’s Formulas yields z∗ 1·z∗ 2=k0 k2 =⇒z∗ 1·z∗ 2<0. 17 Since k1, k2, k3are all real numbers, z∗ 1, z∗ 2must also, so the inequality above must imply that either z∗ 1or z∗ 2is positive, but not both. Let z∗be this unique positive solution to h(z) = 0. Next, we wish to show that 0 < z∗<1. Since R0>1, we know k0<0byLemma 1. Additionally, h(0) = k0=⇒h(0) <0 h(1) = 1 + Th3+Th1+Th3Th1=⇒h(1) >0 since Th1and Th3are greater than 0 when all parameters are biologically possible. Because h(z) is a continuous function with h(0) <0 and h(1) >0, the Intermediate Value Theorem implies that there exists a solution to h(z) = 0 between z= 0 and z= 1. Since z∗is the only positive solution to h(z) = 0, we conclude that 0< z∗<1.(25) Recall that z∗represents the fraction of mosquitoes at the equilibrium point E∗. As such, z∗is biologically possible when R0>1 since 0 < z∗<1. Lastly, we wish to show that y∗ 1, y∗ 2, y∗ 3are also between 0 and 1 when all parameters are biologically possible. Recall that Equation 21 tells shows y∗ i=Thiwiz∗ Thiz∗+ 1 Since wirepresents the fraction of the population with genotype i,wimust satisfy 0 ≤wi≤1 in the real world. As such, the denominator for y∗ iis necessarily greater than the numerator. Therefore, at the nontrivial equilibrium, 0 ≤y∗ i<1, meaning y∗ 1, y∗ 2, and y∗ 3are all biologically possible. As such, we conclude that E∗is biologically possible when all parameters are biologically possible and R0is greater than 1. 2.9 Conditions for Stability of Endemic Equilibrium In the previous section, we proved that the endemic equilibrium E∗exists in biologically possible space only when R0>1. We now present the following theorem on the stability of E∗. Theorem 4. If R0>1and all parameters are biologically possible (i.e. E∗is biologically possible), then the fixed point E∗is stable. Proof. Taking the Jacobian of the fast system (18) and evaluating it at E∗yields J(E∗) =     −(βh1z∗+γ1) 0 0 βh1(w1−y∗ 1) 0−(βh2z∗+γ2) 0 βh2(w2−y∗ 2) 0 0 −(βh3z∗+γ3)βh3(w3−y∗ 3) βv1(1 −z∗)βv2(1 −z∗)βv3(1 −z∗)−(βv1y∗ 1+βv2y∗ 2+βv3y∗ 3+δ)    . 18 Note that J(E∗) = M−Dwhere Mand Dare both non-negative matrices and Dis a diagonal matrix. M=    000βh1(w1−y∗ 1) 000βh2(w2−y∗ 2) 000βh3(w3−y∗ 3) βv1(1 −z∗)βv2(1 −z∗)βv3(1 −z∗) 0    , D=    βh1z∗+γ10 0 0 0βh2z∗+γ20 0 0 0 βh3z∗+γ30 0 0 0 βv1y∗ 1+βv2y∗ 2+βv3y∗ 3+δ    . Recall that Theorem 2 states that the real part of the eigenvalues of J(E∗) are negative if the dominant eigenvalue of MD−1is less than 1. It can be shown that MD−1has a double eigenvalue of 0 and the other two are λ+, λ−=±sA1βv1(1 −z∗) βh1z∗+γ1+A2βv2(1 −z∗) βh2z∗+γ2+A3βv3(1 −z∗) βh3z∗+γ3(26) where Aj=βhj(wj−y∗ j) P3 i=1 [βviy∗ i] + δ(27) Note that λ−is always less than 1 provided that all parameters are within biologically possible bounds. Next, we wish to show that λ+is also less than 1. From Equation 20, recall that 0 = (1 −z∗)(βv1y∗ 1+βv2y∗ 2+βv3y∗ 3)−δz∗. Manipulating this equation to isolate z∗yields z∗=P3 i=1 [βviy∗ i] P3 i=1 [βviy∗ i]−δ. Using the equations of form γjy∗ j=βhjz∗(wj−y∗ j) from Equation 20, we obtain the following expression for z∗. z∗=γjy∗ j βhjz∗(wj−y∗ j), j ∈ {1,2,3} Setting these two expressions for z∗equal and multiplying both sides by the fraction βhj(wj−y∗ j) P3 i=1 [βviy∗ i] yields βhj(wj−y∗ j) P3 i=1 [βviy∗ i]−δ=γjy∗ j P3 i=1 [βviy∗ i], j ∈ {1,2,3} 19 Note that the left side of this expression is equivalent to our definition of Aj,Equation(27). As such, Aj=γjy∗ j P3 i=1 [βviy∗ i], j ∈ {1,2,3}. Furthermore, Ajβvj(1 −z∗) βhjz∗+γj=γjy∗ j P3 i=1 [βviy∗ i]βvj(1 −z∗) βhjz∗+γj, j ∈ {1,2,3} =βvjy∗ j(1 −z∗) P3 i=1 [βviy∗ i]γj βhjz∗+γj. From here, since all parameters are assumed to be biologically possible, βhj +γj> γj. Using this inequality, we find Ajβvj(1 −z∗) βhjz∗+γj<βvjy∗ j(1 −z∗) P3 i=1 [βviy∗ i], j ∈ {1,2,3}.(28) Next, we sum Equation (28) over j= 1,2,3 to obtain the following. Note that the left side is equivalent to (λ+)2by (26). 3 X j=1 Ajβvj(1 −z∗) βhjz∗+γj< 3 X j=1 "βvjy∗ j(1 −z∗) P3 i=1 [βviy∗ i]#j∈ {1,2,3} (λ+)2<1−z∗ Recall from Equation (25) that 0 < Z∗<1. Therefore λ+<1. Since both λ−and λ+are less than 1, Theorem 2 applies, meaning that the real part of the eigenvalues of J(E∗) are negative. As such, we conclude that E∗is stable if R0is greater than 1 and all parameters are biologically possible. 3 Sensitivity Analysis We conduct global sensitivity analysis on our model to determine which parameters are most influential to the model’s behavior. We include biting rate seasonality on our global sensitivity analysis of the fast subsystem, as outlined by the inter-compartmental approach from Renardy et al [21]. There are multiple methods to perform global sensitivity analysis. For nonlinear and monotonic relationships, measures that work the best are Spearman rank correlation coefficient (RCC or Spearman’s ρ), partial rank correlation coefficient (PRCC), and standardized rank regression coefficients (SRRC) [17]. The Sobol method is effective for nonlinear nonmonotonic trends, as do the Fourier amplitude sensitivity test (FAST) and its extended version (eFAST). The derivative based one at a time (OAT) method can be used on any continuous system where it isn’t computationally difficult to calculate the partial derivatives. Since the relationship in the fast subsystem with seasonality is nonlinear and nonmonotonic with computationally possible derivatives we can either utilize the eFAST or the 20 OAT method. We proceed with the OAT method at this time but will conduct eFAST in the future. 3.1 Calculating The Elasticity of R0 We compute the sensitivity indices of the parameters to determine how much each parameter influences R0. The sensitivity index for a parameter pis calculated as ∂R0 ∂p p R0 . The sensitivity indices are as shown below. A contour plot of the values of R0depending on the percentages of the population with XX and XO can be seen in Figure 3. Since the biting rate aand the mosquito death rate δrely on the population of mosquitoes, R0is most drastically influenced by these parameters. The parameters θ3,ϕ3, and γ3all have low elasticities because of the low percentage of the population with absence of the Duffy trait. Figure 3: Sensitivity indices for R0for each parameter in Rio Pardo 3.2 Global Sensitivity on Fast System with Biting Rate Seasonality The elasticity of R0substantiates that the mosquito biting rate is the most influential parameter to R0. We can then account for the seasonality of biting rate by equating ato the sinusoidal model determined by Iyaniwura et al. [15]: a=a0[1 + ϵcos(2πω(t−tc))] . 21 We use derivative-based sensitivity analysis and calculate the normalized sensitivity index from Arriola and Hyman [3] for the equation denoting the total change in the fraction of people getting infected. Thus, we are calculating ∂y ∂p p yfor each parameter pin ywhere y=y1+y2+y3where the other parameters are the baseline values determined by literature. The average sensitivity indices for the Rio Pardo and Macapa DARC distribution are shown in Fig. 4. Macapa has a larger percentage of OO people than Rio Pardo. Therefore, w3has a greater sensitivity index in Macapa but is still smaller than w1and w2due to the reduced transmission rates for OO people. Figure 4: Mean sensitivity index of each parameter. The left graph uses the initial conditions and DARC distribution in Rio Pardo, Brazil while the right graph uses the values from Macapa, Brazil. 4 Parameter Estimation To estimate the transmission rates of malaria between mosquitoes and humans, we conduct parameter estimation on the fast subsystem using incidence data from Brazil and by assuming that Rio Pardo’s Duffy genotype frequency is representative of the wider Brazilian Amazon. We begin by splitting up time series of cases in Brazil from August 2020 to May 2023 from Garcia et al. [9] to the number of people in each DARC genotype infected. 4.1 Relative Risk of Infection Based on Genotype Due to the limitations of available literature, we use our known information to infer the infection rates of P. vivax malaria, based on the Duffy polymorphisms in Rio Pardo. For our model, we attempt to weight the relative risk of being infected with P. vivax malaria (ni) by an individual’s genotype. We presume that if everyone in the population contracts malaria at the same rate, then vi v1 is equal to wi w1 . We used the values of variables provided in the paper by Kano et al. [16] to adjust the risk according to the genotype of each population, resulting in this system of equations: 22            v2 v1 =n2 n1w2 w1 v3 v1 =n3 n1w3 w1 v1+v2+v3=vtotal (29) where niis the relative risk of infection of P. vivax malaria for someone with a genotype i. By calculating v2 v1 and v3 v1 and solving the system of equations, we can determine the proportion of the infected population by genotype. The values of nimust be calculated to determine the relative risk of infection for each genotype. As determined by Kano et al. [16], let naa be the risk of infection due to the FY*A/FY*A genotype; nbb be the risk of infection due to the FY*B/FY*B genotype: n1=FY*A/FY*A FY*A/FY*A + FY*B/FY*B + FY*A/FY*B(naa)+ FY*B/FY*B FY*A/FY*A + FY*B/FY*B + FY*A/FY*B(nbb)+ FY*A/FY*B FY*A/FY*A + FY*B/FY*B + FY*A/FY*B(nab) =160 160 + 93 + 204(0.92) + 93 160 + 93 + 204(1.00) + 204 160 + 93 + 204(1.00) = 0.9719912 (30) We define n2as the weighted average of the relative risk of infection a heterozygous genotype with a silent allele. Let nao be the risk of infection for FY*A/FY*O genotype, and nbo be the risk of the FY*B/FY*O genotype. n2=FY*A/FY*O FY*A/FY*O + FY*B/FY*O(nao) +FY*B/FY*O FY*A/FY*O + FY*B/FY*O(nbo), =139 139 + 73(0.81) + 73 139 + 73(1.26), = 0.9649528. (31) The relative risk n2 n1 of the heterozygous genotype is therefore calculated as follows: n2 n1 =0.9649528 0.9719912 = 0.992758.(32) Thus, from Equation 29 we solve for v2 v1 and v3 v1 : 23 v2 v1 =n2 n1w2 w1v3 v1 =n3 n1w3 w1 = (0.9649528) 0.3072 0.6623= (0.09) 0.0304 0.6623 = 0.4475819 = 0.00413 where v3 v1 represents the homozygous Duffy-negative genotype. 4.2 Fitting the Fast System with Seasonality We begin the parameter estimation by fitting the sinusoidal seasonality of biting rate by fitting the below model from [15] to the incidence cases of Brazil: a=a0h1 + ϵcos 2πω(t−tc)i. Dual annealing was used to minimize the mean absolute percent error (MAPE) between aand the incidence of cases. This was used to fit the parameter ω. We then continue to estimate every parameter except ωwith dual annealing where it minimizes the following value: MAPE(Total incidence) = MAPE(Incidence of XX cases) + MAPE(Incidence of XO cases) + MAPE(Incidence of OO cases). Figure 5: Modeled mosquito biting rate and P. vivax incidence in the Brazilian Amazon 24 Figure 6: Fit of model to the estimated breakdown of P. vivax by Duffy positivity in the Brazilian Amazon In Fig. 5, the biting rate lags the incidence of vivax malaria cases. This is because it takes time for an increased number of bites to result in more infected mosquitoes, which would in turn increase in the number of infected people. In Fig. 6, the lines representing the data look very similar because they differ by a constant factor due to our assumptions when disaggregating the incidence data. Likewise, the model looks very similar because at any given time, the rate of Duffy-positive individuals become infected differs from the rate of infection of Duffy-negative individuals. We assumed that the distribution of people with each DARC genotype was equal to the distribution in Rio Pardo. The results of the parameter estimation can be seen in Table 4. Table 3: Estimated Parameters for Fast System with Biting Rate Seasonality Human Parameters Description Value θ1Probability of infection from infected bite for XX 0.729 θ2Probability of infection from infected bite for XO 0.729 θ3Probability of infection from infected bite for OO 0.060 γ1Recovery Rate for XX 0.0159 γ2Recovery Rate for XO 0.0159 γ3Recovery Rate for OO 0.0159 Mosquito Parameters a0Average Biting Rate 0.0613 25 References [1] Andargie Abate et al. “Differential transmissibility to Anopheles arabiensis of Plasmodium vivax gametocytes in patients with diverse Duffy blood group genotypes”. In: Malaria Journal 22.1 (2023), p. 136. [2] FB Agusto, AB Gumel, and PE Parham. “Qualitative assessment of the role of temperature variations on malaria transmission dynamics”. In: Journal of Biological Systems 23.04 (2015), p. 1550030. [3] Leon Arriola and James M Hyman. “Sensitivity analysis for uncertainty quantification in mathematical models”. In: Mathematical and statistical estimation approaches in epidemiology. Springer, 2009, pp. 195–247. [4] Brazil. en. https://data.who.int/countries/076. Accessed: 2024-6-19. [5] Thomas S Churcher et al. “Probability of transmission of malaria from mosquito to human is regulated by mosquito parasite density in naive and vaccinated hosts”. In: PLoS pathogens 13.1 (2017), e1006108. [6] Zhilan Feng and Carlos Castillo-Chavez. “The influence of infectious diseases on population genetics”. In: Mathematical Biosciences & Engineering 3.3 (2006), pp. 467– 483. [7] Zhilan Feng et al. “Coupling ecology and evolution: malaria and the S-gene across time scales”. In: Mathematical biosciences 189.1 (2004), pp. 1–19. [8] Renzo Galanello and Raffaella Origa. “Beta-thalassemia”. In: Orphanet journal of rare diseases 5 (2010), pp. 1–15. [9] Klauss Kleydmann Sabino Garcia et al. “Is Brazil reaching malaria elimination? A time series analysis of malaria cases from 2011 to 2023”. In: PLOS Global Public Health 4.1 (2024), e0002845. [10] Carlos A Guerra et al. “The international limits and population at risk of Plasmodium vivax transmission in 2009”. In: PLOS neglected tropical diseases 4.8 (2010), e774. [11] Terence J Hadley and Stephen C Peiper. “From malaria to chemokine receptor: the emerging physiologic role of the Duffy blood group antigen”. In: Blood, The Journal of the American Society of Hematology 89.9 (1997), pp. 3077–3091. [12] Gabriela H¨oher, Marilu Fiegenbaum, and Silvana Almeida. “Molecular basis of the Duffy blood group system”. In: Blood Transfusion 16.1 (2018), p. 93. [13] Rosalind E Howes et al. “Global epidemiology of Plasmodium vivax”. In: The American journal of tropical medicine and hygiene 95.6 Suppl (2016), p. 15. [14] Rosalind E Howes et al. “The global distribution of the Duffy blood group”. In: Nature communications 2.1 (2011), p. 266. [15] Sarafa Adewale Iyaniwura et al. “Regional variation and epidemiological insights in malaria underestimation in Cameroon”. In: medRxiv (2023), pp. 1–11. [16] Flora Satiko Kano et al. “Susceptibility to Plasmodium vivax malaria associated with DARC (Duffy antigen) polymorphisms is influenced by the time of exposure to malaria”. In: Scientific reports 8.1 (2018), p. 13851. 32 [17] Simeone Marino et al. “A methodology for performing global uncertainty and sensitivity analysis in systems biology”. In: Journal of theoretical biology 254.1 (2008), pp. 178– 196. [18] GM1 Meny. “The Duffy blood group system: a review”. In: Immunohematology 26.2 (2010), pp. 51–56. [19] Joseli Oliveira-Ferreira et al. “Malaria in Brazil: an overview”. In: Malaria journal 9 (2010), pp. 1–15. [20] Jean Popovici, Camille Roesch, and Virginie Rougeron. “The enigmatic mechanisms by which Plasmodium vivax infects Duffy-negative individuals”. In: PLoS Pathogens 16.2 (2020), e1008258. [21] Marissa Renardy et al. “Global sensitivity analysis of biological multiscale models”. In: Current opinion in biomedical engineering 11 (2019), pp. 109–116. [22] David J Sullivan, Peter Agre, et al. “Human Plasmodium vivax mosquito experimental transmission”. In: The Journal of Clinical Investigation 130.6 (2020), pp. 2800–2802. [23] Priya Venkatesan. “The 2023 WHO World malaria report”. In: The Lancet Microbe 5.3 (2024), e214. [24] Amy Yomiko Vittor et al. “The effect of deforestation on the human-biting rate of Anopheles darlingi, the primary vector of falciparum malaria in the Peruvian Amazon”. In: American Journal of Tropical Medicine and Hygiene 74.1 (2006), pp. 3–11. [25] Liping Wang et al. “Modeling the transmission and control of Zika in Brazil”. In: Scientific reports 7.1 (2017), p. 7721. [26] Michael T White et al. “Immunogenicity of the RTS, S/AS01 malaria vaccine and implications for duration of vaccine efficacy: secondary analysis of data from a phase 3 randomised controlled trial”. In: The Lancet infectious diseases 15.12 (2015), pp. 1450– 1458. [27] Michael T White et al. “Plasmodium vivax and Plasmodium falciparum infection dynamics: re-infections, recrudescences and relapses”. In: Malaria journal 17 (2018), pp. 1–15. 33