A MALE-FEMALE MATHEMATICAL MODEL OF HUMAN PAPILLOMAVIRUS (HPV) IN AFRICAN AMERICAN POPULATION

. We introduce mathematical human papillomavirus (HPV) epidemic models (with and without vaccination) for African American females (AAF) and African American males (AAM) with “ﬁtted” logistic demographics and use these models to study the HPV disease dynamics. The US Census Bureau data of AAF and AAM of 16 years and older from 2000 to 2014 is used to “ﬁt” the logistic demographic models. We compute the basic reproduction number, R 0 , and use it to show that R 0 is less than 1 in the African American (AA) population with or without implementation of HPV vaccination program. Furthermore, we obtain that adopting a HPV vaccination policy in the AAF and AAM populations lower R 0 and the number of HPV infections. Sensitivity analysis is used to illustrate the impact of each model parameter on the basic reproduction number.


1.
Introduction. In the United States of America, human papillomavirus (HPV) is the most common sexually transmitted infection (STI) in males and females [8]. Most sexually active males and females will get at least one type of HPV infection at some point in their lives [5]. In the United States, about 79 million are currently infected with HPV and about 14 million people become newly infected each year [8]. There are more than 150 different types of HPV [7]. Health problems related to HPV include genital warts and cancer. Most people infected with genital HPV do not know they are infected and never develop symptoms or health problems from it. Some people find out they have HPV when they get genital warts. Females may find out they have HPV when they get an abnormal Pap test result during cervical cancer screening. Others may only find out once they have developed more serious problems from HPV, such as cancer [5]. Most HPV infections cause no symptoms and are not clinically significant, but persistent infection can lead to disease or cancer.
Mathematical epidemic models have been used to study HPV infections in various populations. For example, Alsaleh and Gummel [1] in a recent paper, used a deterministic model to assess the impact of vaccination on both high-risk and lowrisk HPV infection types. Ribassin-Majed and Clemencon [14] used a deterministic mathematical model to assess the impact of vaccination on non-cancer causing HPV (6/11) in French males and females. Lee and Tameru [13] used a deterministic model to assess the impact of HPV on cervical cancer in African American females (AAF).
In all these studies, the HPV models have a constant recruitment function for the demographic equations.
In this paper, we use a two-sex HPV model with "fitted" logistic demographics to study the HPV disease dynamics in AAF and African American males (AAM) of 16 years and older. Using US Census Bureau data for AAF and AAM populations, we illustrate that the "fitted" logistic demographic equation captures the African American (AA) population better than the constant recruitment demographic equation. We compute the basic reproduction number, R 0 , and perform sensitivity analysis on R 0 . We obtain that in the AA population R 0 < 1. In addition, we use an extension of the model with vaccination classes to assess the impact of vaccination on the AAF and AAM populations.
The paper is organized as follows: In Section 2, we introduce a demographic equation for AAF (respectively, AAM) and we "fit" it to the US Census Bureau data of AAF (respectively, AAM) of 16 years and older. We introduce, in Section 3, a two-sex African American HPV model. In Section 4, we study disease-free equilibria and compute the basic reproduction number R 0 . In Section 5, we introduce the model with vaccination. In Section 6, we study disease-free equilibria and compute the basic reproduction number R v 0 for HPV model with vaccination. We summarize our results in Section 7.
2. Demographic equations. In [14], Ribassin-Majed et al. used a HPV model with constant recruitment rate in the demographic equation to study HPV disease dynamics in male and female populations of France. In the absence of the HPV disease, the demographic equation of their model is the following ordinary differential equation: where N is total population of males or females and for both males and females µ is the constant per capita mortality rate and Λ is the constant recruitment rate. In [13], Lee and Tameru used Model (1) as the demographic equation for their single sex HPV model in African Americans.
In this paper, we use logistic models that are "fitted" to the 2000 − 2014 US Census Bureau population data of AAF and AAM of 16 years and older (see Tables  1-2) as the demographic equations for AAF and AAM in a HPV model. From Tables 1 and 2, we note that both male and female populations of African Americans of 16 years and older as well as the total populations exhibit increasing trends from 2000 to 2014. In the next section, we use a logistic differential equation model to "capture" the AA population data of Tables 1 and 2.
2.1. Demographic equation for AAF and AAM. In [3], Brauer and Castillo-Chavez used the logistic equation "fitted" to United States Census Bureau data to model the total United States population. We use the same approach to "fit" the solution of the following logistic equation to the AAF population of 16 years and older (see Table 1) and the AAM population of 16 years and older (see Table 2).
where index i = f refers to the female and index i = m refers to the male. N i (t) is the total population of AAF of 16 years and older at time t if i = f , respectively, the total population of AAM of 16 years and older at time t if i = m. The constant µ i is the death rate of population, r i is the intrinsic growth rate of population and K i is the carrying capacity.  Tables 1 and  2 is t = 0. The solutions of (2) are Let R di is the demographic threshold of the population. When the constant intrinsic growth rate, r i , is bigger than the per capita mortality rate of population, µ i , then R di > 1. However, R di < 1 when r i is less than µ i . Consequently, for every N 0 i > 0, when R di > 1 then the nontrivial solution N i (t) approaches the "death" adjusted carrying capacity, and the population persists. However, when R di < 1 then N i (t) → 0 as t → ∞ and the population goes extinct.
2.2. AAF and AAM "fitted" logistic demographic equations. Equation (2) gives the per capita growth rate, dNi/dt Ni , for the total AAF population of 16 years and older when i = f and the total AAM population of 16 years and older when i = m. That is, if N i (t) = 0, then Using the 2002 − 2014 US Census Bureau data of Table 1 (respectively, Table  2), we estimate the values of dN f /dt (respectively, dN m /dt) by symmetric differences. We then graph the per capita growth rate, Y = Fitting the line to the curve gives r f − µ f = 0.021298978 > 0 and r f /K f = 4.02117×10 −10 (respectively, r m −µ m = 0.019877926 > 0 and r m /K m = 1.9980893× 10 −10 ). So, we estimate the intrinsic growth rate of AAF (respectively, AAM) to be r f = 0.028564978 (respectively, r m = 0.028104926) and the carrying capacity of AAF (respectively, AAM) is K f = 71, 036, 484 (respectively, K m = 140, 659, 009). From our estimates of r f and µ f we note that R d f > 1 and K * f = 52, 967, 117 (rounded). Similarly, from our estimates of r m and µ m we note that R dm > 1 and K * m = 99, 484, 673 (rounded). Using our estimates, we express the nontrivial solution (3) of the logistic growth model for AAF of 16 years and older as and for AAM of 16 years and older as N m (t) = 99, 484, 673 where year 2001 is taken as t = 0. The plot of the data of AAF of 16 years or older and solution (5) in Figure 1 show that our "fitted" model captures the AAF data of Table 1. Similarly, the plot of the data of AAM of 16 years or older and solution (6) in Figure 2 show that our "fitted" model captures the AAM data of Table 2. 2.3. AAF and AAM "fitted" logistic equation versus constant recruitment model. When the population in equation (1) consists only of AAF (i = f ) or of AAM (i = f ) of 16 years and older, then the solution is Using the initial condition N 0 f = 14, 041, 520 [12] (respectively, N 0 m = 12, 124, 810 [12]), the parameters Λ f = 16, 821, 072 (52% of AA population are female [2,13]) and µ f = 0.007266 [6] (respectively, Λ m = 15, 527, 714 (48% of AA population are male [2, 13]) and µ m = 0.008227 [6]), solution (7) becomes N f (t) = 2, 315, 038, 811 − 2, 300, 997, 291 e −0.007266t , and N m (t) = 1, 887, 409, 019 − 1, 875, 284, 209 e −0.008227t . (9) In Figure 1, we compare the "fitted" solution (8) of Model (1), and our "fitted" solution (5) of Model (2), to the US Census Bureau data in Table 1. Figure 1 shows that Model (2), the "fitted" logistic model, captures better the 2000 − 2014 US Census Bureau data of the AAF population of 16 years and older than Model (1), the constant recruitment demographic equation used in the France study [14] and in the study of Lee and Tameru [13]. Similarly, Figure 2 shows that, as in the female population, the "fitted" solution (9) of Model (1) over estimates the AAM US Census Bureau data while our "fitted" solution (5) of Model (2) captures it.

3.
A two-sex African American HPV model. To study the HPV dynamics in male and female African American populations of 16 years and older, we assume that the total AAF population (respectively, total AAM population) of 16 years and older is governed by Model (2) with i = f (respectively, Model (2) with i = m). As in [14], we divide the population into four compartments: S f = susceptible AAF population, S m = susceptible AAM population, I f = HPV infected AAF population, I m = HPV infected AAM population. It is known that females and males recover from HPV at about the same rate [11]. We now introduce the following mathematical model that uses our "fitted" logistic demographics for the total population.

344
NAJAT ZIYADI where N f = S f + I f > 0, N m = S m + I m > 0 and the model parameters with their values are listed in Table 3. Infection rate for AAM population [1] Notice that since adding the S f −equation to the I f −equation (respectively, adding the S m −equation to the I m −equation) gives the AAF demographic equation (2) with i = f (respectively, AAM demographic equation (2) with i = m). Consequently, in Model (10), the total AAF and AAM populations, governed by our "fitted" logistic equations (5) and (6) are bounded. We will study Model (10) with the parameter values listed in Table 3 and with the initial conditions listed in Table 4. Notice that in Table 4 Table 1 (10), we obtain that Similarly, if I 0 m = 0, then from the I m −equation of Model (10), we obtain that If S 0 f = 0, then from the S f −equation of Model (10), we obtain that (10), we obtain that Hence, S f and I f are bounded. Similarly, S m and I m are bounded.
Then Model (10) is equivalent to In AAF population, R d f > 1 implies that for t > 0, N f (t) > 0 and lim t→∞ N f (t) = K * f when N f (0) > 0. Similarly, in AAM population, R dm > 1 implies that for t > 0, N m (t) > 0 and lim t→∞ N m (t) = K * m when N m (0) > 0. Hence, N f (t) + N m (t) > 0 in Model (10) and F is C 1 . Consequently, with our parameters from Table 3 and  initial conditions from Table 4, Model (10) has a unique nonnegative solution for t ≥ 0.
4. Disease-free equilibria, stability and basic reproduction number (R 0 ). Unlike in [14], it is possible for Model (10) to exhibit up to four disease-free equilibrium points (DFEs), where R d f > 1 and R dm > 1.
Similarly, to determine the stability of P m0 , we compute the Jacobian matrix at P m0 .
To determine the stability of P f m , we recall that in this case, for N 0 Consequently, Model (10) reduces to the following "limiting" system: Using the next generation matrix method [17] we obtain the following two matrices Then, using the Jacobian matrices of F and V evaluated at the DFE of (11), Q f m = (0, 0, K * f , K * m ), we obtain the following matrices, Hence, By the next generation matrix method [17], the reproduction number for Model (10), R 0 , is the spectral radius of the next generation matrix, ρ(F V −1 ), and . R 0f (respectively, R 0m ) is directly proportional to the product of the infection rate of AAF population and the proportion of AAF population at the death-adjusted equilibrium point (respectively, the infection rate of AAM population and the proportion of AAM population at the death-adjusted equilibrium point) and inversely proportional to the sum of clearance rate and death rate for AAF population (respectively, the sum of clearance rate and death rate for AAM population).
When r f > µ f and r m > µ m , then R 0 < 1 implies the DFE, P f m , is locally asymptotically stable, whereas R 0 > 1 implies P f m is unstable [17]. Using the parameter values in Table 3, it is easy to see that R 0f = 0.1915 < 1 and R 0m = 0.2874 < 1. Hence, R 0 = 0.2346 < 1. That is, in our Model (10) the DFE P f m , is locally asymptotically stable.
It is interesting to note that in [13], Lee and Tameru, obtained R 0 = 0.519798 < 1 with a one sex HPV model under a constant recruitment. Since their model over estimates the AA population, their R 0 is bigger than 0.2346. However, both models predict that R 0 < 1 and HPV infection cannot get started in a fully susceptible AA population [13].
To find an effective mitigation strategy that seeks to reduce HPV infection in AA population within the shortest time possible, in the next section, we use sensitivity analysis to study the impact of each model parameter on R 0 . 4.1. Sensitivity analysis of R 0 . Sensitivity indices are used to measure the relative change in a state variable when a parameter changes. Typically, the normalized forward sensitivity index of a variable to a parameter is defined as the ratio of the relative change in the variable to the relative change in the parameter. When the variable is a differential function of the parameter, the sensitive index may be alternatively defined using partial derivatives [4,10,18,19]. Definition 4.1 ([4,10,18,19]). The normalized forward sensitivity index of a variable, u, that depends differentiably on a parameter, q, is defined as: We use Definition 4.1 to derive the sensitivity indices of the basic reproduction number R 0 and we evaluate them using parameter values of Table 3.
Increasing (respectively, decreasing) the clearance rate, δ, by 1% will decrease (respectively, increase) the value of R 0 by about 0.99%. Increasing (respectively, decreasing) the infection rate of the AAF population, σ f , by 1% will increase (respectively, decrease) the value of R 0 by about 0.5%. Increasing (respectively, decreasing) the infection rate of the AAM population, σ m , by 1% will increase (respectively, decrease) the value of R 0 by about 0.5%.

African American male and female HPV model simulations.
To illustrate the impact of HPV on AAF and AAM populations of 16 years and older, we simulate Model (10) with the parameter values listed in Table 3 and the initial conditions listed in Table 4. Table 5. Normalized sensitivity indices and order of importance of R 0 to the nine parameters in Table (3).  Table (3). The most sensitive parameters for R 0 are the clearance rate, δ, the infection rate of the AAF population, σ f , and the infection rate of the AAM population, σ m . While the least sensitive parameters are the intrinsic growth rate for AAF population, r f , the intrinsic growth rate for AAM population, r m , the death rate of AAF population, µ f , and the death rate for AAM population, µ m .
Simulations of our HPV Model (10) are performed using Matlab software, and are illustrated in Figure 4. Figure 4 (a) shows that susceptible population of AAF of 16 years and older, S f , increases over time. Similarly, Figure 4 (b) shows that susceptible population of AAM of 16 years and older, S m , increases over time. In the total AA population, R 0 < 1 and as expected, Figure 4 (c) shows that the AAF HPV infected population of 16 years and older, I f , decreases monotonically with time and Figure 4 (d) shows that the AAM HPV infected population of 16 years and older, I m , decreases monotonically with time.
To protect against HPV infections, HPV vaccines are available for males and females. Gardasil and Cervarix are two HPV vaccines that have market approval in many countries. Next, we introduce an extension of Model (10) with vaccinated male and female classes. We will use the extended model to study the impact of vaccination on Figure 4.

5.
A two-sex African American HPV model with vaccination. To introduce HPV vaccination in AAF and AAM populations of Model (10), we let p f (respectively, p m ) denote the proportion of HPV vaccinated females (respectively, males). For both males and females, we let τ denote the success rate of the vaccine. In AA population, τ = 90% [14]. As in [14], we divide the AA population into eight compartments. S f = non vaccinated susceptible AAF population, S v f = vaccinated susceptible AAF population, S m = non vaccinated susceptible AAM population, S v m = vaccinated susceptible AAM population, I f = non vaccinated HPV infected AAF population, I v f = vaccinated HPV infected AAF population, I m = non vaccinated HPV infected AAM population, I v m = vaccinated HPV infected AAM population. Then Model (10) with vaccination in both AAF and AAM of 16 years and older becomes the following model.   Table 3. We will study Model (12) with parameter values listed in Table 3 and with the initial conditions listed in Table 6, where p f = 39%, p m = 20.4% and τ = 90%.
Note that in Table 6, Table 1 Table 2) where t = 0 is year 2001.
Proceeding exactly as in Theorem 3.1, we obtain the following result. From Model (12), the demographic equations for the female and male total populations are respectively the following: The equilibrium points of Model (13) are (N f , N m ) = (0, 0), (0, K * m ), (K * f , 0) and (K * f , K * m ). Since R d f > 1 and R dm > 1, (N f , N m ) = (0, 0), (0, K * m ) and (K * f , 0) are unstable and (K * f , K * m ) is asymptotically stable. As in Model (10), to state the "limiting" system of Model (12), we replace N f by K * f and N m by K * m . Since the "limiting" system for Model (12) is the following system of equations.
By the next generation matrix method [17], the reproduction number for Model (12), R v 0 , is the spectral radius of the next generation matrix, ρ(F V −1 ), and and .
Thus, adopting a HPV vaccination program decreases the basic reproduction number, R 0 , in AA population.
Using the parameter values of Table 3, τ = 90% [14], p f = 39% [9] and p m = 20.4% [9], we obtain that R v 0f = 0.1243 and R v 0m = 0.2346. Hence, R v 0 = 0.1708. In the next section, we use sensitivity analysis to illustrate the impact of model parameters on R v 0 . 6.1. Sensitivity analysis of R v 0 . We use Definition 4.1 to derive the sensitivity indices of the basic reproduction number R v 0 and we evaluate them using, τ = 90%, p f = 39%, p m = 20.4% and parameter values of Table (3). Table 7. Normalized sensitivity indices and order of importance of R v 0 to model parameters. are the clearance rate, δ, the infection rate of the AAF population, σ f , and the infection rate of the AAM population, σ m , the success rate of HPV vaccine, τ , and the proportion of HPV vaccinated females, p f . While the least sensitive parameters are the intrinsic growth rate for AAF population, r f , the intrinsic growth rate for AAM population, r m , the death rate of AAF population, µ f , and the death rate for AAM population, µ m . AAM of 16 years and older, S m (respectively, vaccinated AAM of 16 years and older S v m ), increases over time. In the total population of African American population, R v 0 < 1 and as expected, Figures 6 (e-f) show that the AAF HPV infected population of 16 years and older, I f (respectively, vaccinated AAF of 16 years and older I v f ), decreases monotonically with time and Figures 6 (g-h) show that the AAM HPV infected population of 16 years and older, I m (respectively, vaccinated AAM of 16 years and older I v m ), decreases monotonically with time. To study the impact of the presence of the vaccinated class on the results of Figure 4, we simulate Model (12) using initial conditions in Table 8. For these simulations of Model (12), we keep all the parameter values at their current values in Figure 4, where τ = 90%, while we vary p f and p m .
Note that in Figure 6 and Table 6, S f (0) + S v f (0) (respectively, I f (0) + I v f (0)) is the same as S f (0) (respectivrly, I f (0)) of Table 4 and Figure 4. Similarly, in Figure  6 and Table 6, S m (0) + S v m (0) (respectively, I m (0) + I v m (0)) is the same as S m (0) (respectivrly, I m (0)) of Table 4 and  6.3. Impact of HPV vaccination. In AA population, R v 0 < R 0 < 1. Consequently, with and without vaccination, the susceptible population increases while the infective population decreases over time. In Figures 7 and 8, we illustrate that in both AAF and AAM populations, the increase in the susceptible populations is higher when a vaccination policy is adopted than when the population is not being vaccinated.