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 low-risk 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, $\mathcal{R}_0$ , and perform sensitivity analysis on $\mathcal{R}_0$. We obtain that in the AA population $\mathcal{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 $\mathcal{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 $\mathcal{R}_0^v$ 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:

$\frac{dN}{dt}=\Lambda-\mu N\label{eq:l},$

(1)

where $N$ is total population of males or females and for both males and females $\mu$ is the constant per capita mortality rate and $\Lambda$ 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).

$\frac{dN_i}{dt}=(r_i-\mu_i)N_i(1-\frac{N_i}{K_i(r_i-\mu_i)/r_i}), t\geq0,\label{eq:a}$

(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 $\mu_i$ is the death rate of population, $r_i$ is the intrinsic growth rate of population and $K_i$ is the carrying capacity.

$N_i^0$ is the size of the population at time $t=0$, where year 2001 in Tables 1 and 2 is $t=0$. The solutions of (2) are

$N_i(t)=0 \;\;\text{and} \;\;N_i(t)=\frac{K_i(r_i-\mu_i)/r_i}{1+(\frac{K_i(r_i-\mu_i)/r_i}{N_i^0}-1)e^{-(r_i-\mu_i)t}}.\label{eq:b}$

(3)

Let

$R_{d_i}=\frac{r_i}{\mu_i}.$

$R_{d_i}$ 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, $\mu_i$, then $R_{d_i}>1.$ However, $R_{d_i}<1$ when $r_i$ is less than $\mu_i$. Consequently, for every $N^0_i>0$, when $R_{d_i}>1$ then the nontrivial solution $N_i(t)$ approaches the "death" adjusted carrying capacity,

$K_i^*=\frac{K_i(r_i-\mu_i)}{r_i}, \;\;\text{as}\;\; t\rightarrow \infty,$

and the population persists. However, when $R_{d_i}<1$ then $N_i(t)\rightarrow 0$ as $t\rightarrow \infty$ and the population goes extinct.

---

2.2. AAF and AAM "fitted" logistic demographic equations

Equation (2) gives the per capita growth rate, $\frac{dN_i/dt}{N_i}$, 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)\neq 0$, then

$ \frac{dN_i/dt}{N_i}=(r_i-\mu_i)(1-\frac{N_i}{K_i(r_i-\mu_i)/r_i}). $

(4)

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=\frac{dN_f/dt}{N_f}$ (respectively, $Y=\frac{dN_m/dt}{N_m}$) versus the AAF (respectively, AAM) population. Then the line with equation $Y=0.021298978-4.02117\times10^{-10}N_f$ (respectively, $Y= 0.019877926-1.9980893\times10^{-10}N_m$) "fits" the resulting curve. Therefore, the AAF (respectively, AAM) population of 16 years and older exhibits a logistic behavior (see Figures 1 and 2, respectively).

Fitting the line to the curve gives $r_f-\mu_f=0.021298978>0$ and $r_f/K_f=4.02117\times10^{-10}$ (respectively, $r_m-\mu_m=0.019877926>0$ and $r_m/K_m=1.9980893\times10^{-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 $\mu_f$ we note that $R_{d_f}>1$ and $K_f^*=52,967,117$ (rounded). Similarly, from our estimates of $r_m$ and $\mu_m$ we note that $R_{d_m}>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

$N_f(t) = \frac{52,967,117}{1+2.772176873e^{-0.021298978t}},\label{eq:c}$

(5)

and for AAM of 16 years and older as

$N_m(t) = \frac{99,484,673}{1+7.205050059e^{-0.019877926t}},\label{eq:c-m}$

(6)

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

$ N_i(t)=\frac{\Lambda_i}{\mu_i}+(N_i^0-\frac{\Lambda_i}{\mu_i})e^{-\mu_i t}.\label{eq:e} $

(7)

Using the initial condition $N_f^0=14,041,520$ [12] (respectively, $N_m^0=12,124,810$ [12]), the parameters $\Lambda_f=16,821,072$ ($52\%$ of AA population are female [2,13]) and $\mu_f=0.007266$ [6] (respectively, $\Lambda_m=15,527,714$ ($48\%$ of AA population are male [2,13]) and $\mu_m=0.008227$ [6]), solution (7) becomes

$ N_f(t)=2,315,038,811-2,300,997,291\;e^{-0.007266 t},\label{eq:o} $

(8)

and

$ N_m(t)=1,887,409,019-1,875,284,209\;e^{-0.008227 t}.\label{eq:o-m} $

(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.

$ \left\{ $$\begin{aligned}\label{eq:model-HPV} \frac{dS_f}{dt} = r_fN_f(1-\frac{N_f}{K_f}) - \sigma_fS_f\frac{I_m}{N_f+N_m} + \delta I_f - \mu_f S_f,\\ \frac{dS_m}{dt} = r_mN_m(1-\frac{N_m}{K_m}) - \sigma_mS_m\frac{I_f}{N_f+N_m} + \delta I_m - \mu_m S_m,\\ \frac{dI_f}{dt} = \sigma_fS_f\frac{I_m}{N_f+N_m} - (\delta+\mu_f) I_f,\\ \frac{dI_m}{dt} = \sigma_mS_m\frac{I_f}{N_f+N_m} - (\delta+\mu_m) I_m, \end{aligned}$$ \right. $

(10)

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.

Notice that since

$N_f(t)=S_f(t)+I_f(t)\;\;\text{and}\;\;N_m(t)=S_m(t)+I_m(t),$

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, $S_f(0)+I_f(0)=N_f(0)$ of Table 1 (respectively, $S_m(0)+I_m(0)=N_m(0)$ of Table 2), where $t=0$ is year 2001.

---

3.1. Boundedness of Orbits

In this section, we show that Model (10) is well-posed. In particular, we obtain that all orbits are nonnegative and there is no population explosion in Model (10).

Theorem 3.1. All solutions of Model (10) are nonnegative and bounded.

Proof. Consider the following nonnegative initial conditions $S_f(0)=S_f^0\geq0$, $S_m(0)=S_m^0\geq0$, $I_f(0)=I_f^0\geq0$ and $I_m(0)=I_m^0\geq0$.

If $I_f^0=0$, then from the $I_f-$equation of Model (10), we obtain that $\frac{dI_f}{dt}|_{t=0}=\sigma_fS_f^0\frac{I_m^0}{N_f^0+N_m^0}\geq 0.$ Hence, $I_f(t)\geq0$ for all $t\geq0$ [15].

Similarly, if $I_m^0=0$, then from the $I_m-$equation of Model (10), we obtain that $\frac{dI_m}{dt}|_{t=0}=\sigma_mS_m^0\frac{I_f^0}{N_f^0+N_m^0}\geq 0.$ Hence, $I_m(t)\geq0$ for all $t\geq0$ [15].

If $S_f^0=0$, then from the $S_f-$equation of Model (10), we obtain that $\frac{dS_f}{dt}|_{t=0} = r_fI^0_f(1-\frac{I_f^0}{K_f})+ \delta I_f^0$. $I_f^0\leq K_f$. Hence, $S_f(t)\geq0$ for all $t\geq0$.

Similarly, if $S_m^0=0$, then from the $S_m-$equationn of Model (10), we obtain that $\frac{dS_m}{dt}|_{t=0} = r_mI^0_m(1-\frac{I_m^0}{K_m})+ \delta I_m^0$. $I_m^0\leq K_m$. Hence, $S_m(t)\geq0$ for all $t\geq0$.

Recall that $N_f(t)$ is bounded, $S_f(t)+I_f(t)=N_f(t)$, $S_f(t)\geq0$ and $I_f(t)\geq0$. Hence, $S_f$ and $I_f$ are bounded. Similarly, $S_m$ and $I_m$ are bounded.

Let $X= $$\begin{bmatrix}S_f\\S_m\\I_f\\I_m\end{bmatrix}$$ $ and $F(X)= $$\begin{bmatrix} r_fN_f(1-\frac{N_f}{K_f}) - \sigma_fS_f\frac{I_m}{N_f+N_m} + \delta I_f - \mu_f S_f\\ r_mN_m(1-\frac{N_m}{K_m}) - \sigma_mS_m\frac{I_f}{N_f+N_m} + \delta I_m - \mu_m S_m\\ \sigma_fS_f\frac{I_m}{N_f+N_m} - (\delta+\mu_f) I_f\\\sigma_mS_m\frac{I_f}{N_f+N_m} - (\delta+\mu_m) I_m\end{bmatrix}$$ .$

Then Model (10) is equivalent to

$\frac{dX}{dt}=F(X).$

In AAF population, $R_{d_f}>1$ implies that for $t>0$, $N_f(t)>0$ and $\lim_{t\rightarrow\infty}N_f(t)=K^*_f$ when $N_f(0)>0$. Similarly, in AAM population, $R_{d_m}>1$ implies that for $t>0$, $N_m(t)>0$ and $\lim_{t\rightarrow\infty}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\geq 0$.

---

4. Disease-free equilibria, stability and basic reproduction number ($\mathcal{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_{d_m}>1$.

Since $R_{d_f}>1$ and $R_{d_m}>1$, then Model (10) has four DFEs, $P_{00}=(0,0,0,0)$, $P_{f0}=(K_f^*,0,0,0)$, $P_{m0}=(0,K_m^*,0,0)$, and a DFE with both susceptible AAF and AAM with no HPV infected individuals, $P_{fm}=(K_f^*,K_m^*,0,0).$

Since $R_{d_f}>1$ and $R_{d_m}>1$, $P_{00}$ is unstable.

To determine the stability of $P_{f0}$, we compute the Jacobian matrix at $P_{f0}$.

$J|_{P_{f0}}= $$\begin{bmatrix} \mu_f -r_f&0& 2\mu_f-r_f+\delta &-\sigma_f\\ 0 & r_m-\mu_m& 0 & 2\mu_m-r_m+\delta\\ 0 & 0& -(\delta+\mu_f)& \sigma_f \\ 0 &0 & 0& -(\delta+\mu_m)\\ \end{bmatrix}$$ .$

$J|_{P_{f0}}$ has the following four eigenvalues: $\lambda_1=\mu_f- r_f$, $\lambda_2=r_m-\mu_m$, $\lambda_3=-(\sigma+\mu_f)$ and $\lambda_4=-(\sigma+\mu_m)$. All model parameters are positive so $\lambda_3<0$ and $\lambda_4<0$. Hence, the stability of $P_{f0}$ is determined by the sign of $\mu_f- r_f$ and $r_m-\mu_m$. $R_{d_m}>1$ implies $r_m>\mu_m>0$ and $\lambda_2>0$. Hence, $P_{f0}$ is unstable.

Similarly, to determine the stability of $P_{m0}$, we compute the Jacobian matrix at $P_{m0}$.

$J|_{P_{m0}}= $$\begin{bmatrix} r_f-\mu_f &0&r_f+\delta &0\\ 0 & \mu_m -r_m& -\sigma_m&2\mu_m-r_m+\delta\\ 0 & 0& -(\delta+\mu_f)& 0 \\ 0 &0 & \sigma_m& -(\delta+\mu_m)\\ \end{bmatrix}$$ .$

$J|_{P_{m0}}$ has the following five eigenvalues: $\lambda_1= r_f-\mu_f$, $\lambda_2=r_m-\mu_m$, $\lambda_3=-(\sigma+\mu_f)$ and $\lambda_4=-(\sigma+\mu_m)$. Hence, the stability of $P_{m0}$ is determined by the sign of $r_f-\mu_f $ and $r_m-\mu_m$. $R_{d_f}>1$ implies $r_f>\mu_f>0$ and $\lambda_1>0$. Hence, $P_{m0}$ is unstable.

To determine the stability of $P_{fm}$, we recall that in this case, for $N_f^0>0$ and $N_m^0>0$, $\underset{t \to \infty}{\lim}N_f(t)=K_f^*$ and $\underset{t \to \infty}{\lim}N_m(t)=K_m^*$. Consequently, Model (10) reduces to the following "limiting" system:

$ \left\{ $$\begin{aligned}\label{eq:limit2} &\frac{dI_f}{dt} =\sigma_f(K_f^*-I_f)\frac{I_m}{K_f^*+K_m^*} -(\delta+\mu_f) I_f,\\ &\frac{dI_m}{dt} =\sigma_m(K_m^*-I_m)\frac{I_f}{K_f^*+K_m^*} -(\delta+\mu_m) I_m,\\ &\frac{dS_f}{dt} = r_fK_f^*(1-\frac{K_f^*}{K_f}) - \sigma_fS_f\frac{I_m}{K_f^*+K_m^*} + \delta I_f - \mu_f S_f,\\ &\frac{dS_m}{dt} = r_mK_m^*(1-\frac{K_m^*}{K_m}) - \sigma_mS_m\frac{I_f}{K_f^*+K_m^*} + \delta I_m - \mu_m S_m. \end{aligned}$$ \right. $

(11)

Using the next generation matrix method [17] we obtain the following two matrices

$ \mathcal{F}= $$\begin{bmatrix} \sigma_f(K_f^*-I_f)\frac{I_m}{K_f^*+K_m^*}\\ \sigma_m(K_m^*-I_m)\frac{I_f}{K_f^*+K_m^*}\\ 0\\ 0 \end{bmatrix}$$ $

$\text{and}\;\;\mathcal{V}= $$\begin{bmatrix} (\delta+\mu_f) I_f\\ (\delta+\mu_m) I_m\\ -r_fK_f^*(1-\frac{K_f^*}{K_f}) + \sigma_fS_f\frac{I_m}{K_f^*+K_m^*} - \delta I_f + \mu_f S_f\\ -r_mK_m^*(1-\frac{K_m^*}{K_m}) + \sigma_mS_m\frac{I_f}{K_f^*+K_m^*} - \delta I_m + \mu_m S_m\end{bmatrix}$$ .$

Then, using the Jacobian matrices of $\mathcal{F}$ and $\mathcal{V}$ evaluated at the DFE of (11), $Q_{fm}=(0,0,K_f^*,K_m^*)$, we obtain the following matrices,

$D\mathcal{F}(Q_{fm})=\left[ $$\begin{array}{c|c} F&0\\ \hline 0&0 \end{array}$$ \right]\;\;\text{and}\;\;D\mathcal{V}(Q_{fm})=\left[ $$\begin{array}{c|c} V&0\\ \hline W&U \end{array}$$ \right],$

where $F= $$\begin{bmatrix} 0&\frac{\sigma_f K_f^*}{K_f^*+K_m^*}\\ \frac{\sigma_m K_m^*}{K_f^*+K_m^*}&0 \end{bmatrix}$$ , V= $$\begin{bmatrix} \delta+\mu_f&0\\ 0&\delta+\mu_m \end{bmatrix}$$ , W= $$\begin{bmatrix} -\delta&\frac{\sigma_f K_f^*}{K_f^*+K_m^*}\\ \frac{\sigma_m K_m^*}{K_f^*+K_m^*}&-\delta \end{bmatrix}$$ $ $\text{and}\;\;U= $$\begin{bmatrix} \mu_f&0\\ 0&\mu_m \end{bmatrix}$$ .$

Hence,

$FV^{-1}=\frac{1}{K_f^*+K_m^*} $$\begin{bmatrix} 0&\frac{\sigma_f K_f^*}{\delta+\mu_m}\\ \frac{\sigma_m K_m^*}{\delta+\mu_f}&0 \end{bmatrix}$$ .$

By the next generation matrix method [17], the reproduction number for Model (10), $\mathcal{R}_0$, is the spectral radius of the next generation matrix, $\rho(FV^{-1})$, and

$\mathcal{R}_0=\rho(FV^{-1})=\sqrt{\mathcal{R}_{0f}\mathcal{R}_{0m}},\text{where }$

$\mathcal{R}_{0f}=\frac{\sigma_f K_f^*}{(\delta+\mu_f)(K_f^*+K_m^*)}\;\;\text{and}\;\;\mathcal{R}_{0m}=\frac{\sigma_m K_m^*}{(\delta+\mu_m)(K_f^*+K_m^*)}.$

$\mathcal{R}_{0f}$ (respectively, $\mathcal{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>\mu_f$ and $r_m>\mu_m$, then $\mathcal{R}_0<1$ implies the DFE, $P_{fm}$, is locally asymptotically stable, whereas $\mathcal{R}_0>1$ implies $P_{fm}$ is unstable [17]. Using the parameter values in Table 3, it is easy to see that $\mathcal{R}_{0f}=0.1915<1$ and $\mathcal{R}_{0m}=0.2874<1$. Hence, $\mathcal{R}_0=0.2346<1$. That is, in our Model (10) the DFE $P_{fm}$, is locally asymptotically stable.

It is interesting to note that in [13], Lee and Tameru, obtained $\mathcal{R}_{0}=0.519798<1$ with a one sex HPV model under a constant recruitment. Since their model over estimates the AA population, their $\mathcal{R}_{0}$ is bigger than $0.2346$. However, both models predict that $\mathcal{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 $\mathcal{R}_{0}$.

---

4.1. Sensitivity analysis of $\mathcal{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:

$\Upsilon^u_q:=\frac{\partial u}{\partial q}\times\frac{q}{u}.$

We use Definition 4.1 to derive the sensitivity indices of the basic reproduction number $\mathcal{R}_{0}$ and we evaluate them using parameter values of Table 3.

Increasing (respectively, decreasing) the clearance rate, $\delta$, by $1\%$ will decrease (respectively, increase) the value of $\mathcal{R}_{0}$ by about $0.99\%$. Increasing (respectively, decreasing) the infection rate of the AAF population, $\sigma_f$, by $1\%$ will increase (respectively, decrease) the value of $\mathcal{R}_{0}$ by about $0.5\%$. Increasing (respectively, decreasing) the infection rate of the AAM population, $\sigma_m$, by $1\%$ will increase (respectively, decrease) the value of $\mathcal{R}_{0}$ by about $0.5\%$.

---

4.2. 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.

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, $\mathcal{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 $\tau$ denote the success rate of the vaccine. In AA population, $\tau=90\%$ [14].

As in [14], we divide the AA population into eight compartments. $S_f =$ non vaccinated susceptible AAF population, $S_f^v =$ vaccinated susceptible AAF population, $S_m =$ non vaccinated susceptible AAM population, $S_m^v =$ vaccinated susceptible AAM population, $I_f =$ non vaccinated HPV infected AAF population, $I_f^v =$ vaccinated HPV infected AAF population, $I_m =$ non vaccinated HPV infected AAM population, $I_m^v =$ vaccinated HPV infected AAM population. Then Model (10) with vaccination in both AAF and AAM of 16 years and older becomes the following model.

$ \left\{ $$\begin{aligned}\label{eq:model-HPV-V} \frac{dS_f}{dt} = r_f(1-p_f)N_f(1-\frac{N_f}{K_f}) - \sigma_fS_f\frac{I_m+I_m^v}{N_f+N_m} + \delta I_f - \mu_f S_f,\\ \frac{dS_f^v}{dt} = p_fr_fN_f(1-\frac{N_f}{K_f}) - (1-\tau)\sigma_fS_f^v\frac{I_m+I_m^v}{N_f+N_m} + \delta I_f^v - \mu_f S_f^v,\\ \frac{dS_m}{dt} = r_m(1-p_m)N_m(1-\frac{N_m}{K_m}) - \sigma_mS_m\frac{I_f+I_f^v}{N_f+N_m} + \delta I_m - \mu_m S_m,\\ \frac{dS_m^v}{dt} = p_mr_mN_m(1-\frac{N_m}{K_m}) - (1-\tau)\sigma_mS_m^v\frac{I_f+I_f^v}{N_f+N_m} + \delta I_m^v - \mu_m S_m^v,\\ \frac{dI_f}{dt} = \sigma_fS_f\frac{I_m+I_m^v}{N_f+N_m} - (\delta+\mu_f) I_f,\\ \frac{dI_f^v}{dt} =(1-\tau) \sigma_fS_f^v\frac{I_m+I_m^v}{N_f+N_m} - (\delta+\mu_f) I_f^v,\\ \frac{dI_m}{dt} = \sigma_mS_m\frac{I_f+I_f^v}{N_f+N_m} - (\delta+\mu_m) I_m,\\ \frac{dI_m^v}{dt} =(1-\tau) \sigma_mS_m^v\frac{I_f+I_f^v}{N_f+N_m} - (\delta+\mu_m) I_m^v,\\ \end{aligned}$$ \right. $

(12)

where $N_f=S_f+S_f^v+I_f+I_f^v>0$ and $N_m=S_m+S_m^v+I_m+I_m^v>0$ and the model parameters and their values are listed in 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 $\tau=90\%$.

Note that in Table 6, $S_f(0)+S_f^v(0)+I_f(0)+I_f^v(0)=N_f(0)$ of Table 1 (respectively, $S_m(0)+S_m^v(0)+I_m(0)+I_m^v(0)=N_m(0)$ of Table 2) where $t=0$ is year 2001.

Proceeding exactly as in Theorem 3.1, we obtain the following result.

Theorem 5.1. All solutions of Model (12) are nonnegative and bounded.

---

6. Disease-free equilibrium, stability and basic reproduction number $\mathcal{R}_0^v$

Notice that when all vaccinated classes are missing $(S_f^v=S_m^v=I_f^v=I_m^v=0)$ and $p_f=p_m=\tau=0$, then Model (12) reduces to Model (10).

From Model (12), the demographic equations for the female and male total populations are respectively the following:

$ \left\{ $$\begin{aligned}\label{eq:model-HPV-V-demogr} &\frac{dN_f}{dt} = r_f N_f(1-\frac{N_f}{K_f}) - \mu_f N_f,\\ &\frac{dN_m}{dt} = r_m N_m1-\frac{N_m}{K_m}) - \mu_m N_m.\\ \end{aligned}$$ \right. $

(13)

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_{d_m}>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 $N_f=S_f+S_f^v+I_f+I_f^v=K_f^*$ and $N_m=S_m+S_m^v+I_m+I_m^v=K_m^*$, the "limiting" system for Model (12) is the following system of equations.

$ \left\{ $$\begin{aligned}\label{eq:model-HPV-V-limit} \frac{dS_f}{dt} = r_f(1-p_f)K_f^*(1-\frac{K_f^*}{K_f}) - \sigma_fS_f\frac{I_m+I_m^v}{K_f^*+K_m^*} + \delta I_f - \mu_f S_f,\\ \frac{dS_m}{dt} = r_m(1-p_m)K_m^*(1-\frac{K_m^*}{K_m}) - \sigma_mS_m\frac{I_f+I_f^v}{K_f^*+K_m^*} + \delta I_m - \mu_m S_m,\\ \frac{dI_f}{dt} = \sigma_fS_f\frac{I_m+I_m^v}{K_f^*+K_m^*} - (\delta+\mu_f) I_f,\\ \frac{dI_f^v}{dt} =(1-\tau) \sigma_f(K_f^*-S_f-I_f-I_f^v)\frac{I_m+I_m^v}{K_f^*+K_m^*} - (\delta+\mu_f) I_f^v,\\ \frac{dI_m}{dt} = \sigma_mS_m\frac{I_f+I_f^v}{K_f^*+K_m^*} - (\delta+\mu_m) I_m,\\ \frac{dI_m^v}{dt} =(1-\tau) \sigma_m(K_m^*-S_m-I_m-I_m^v)\frac{I_f+I_f^v}{K_f^*+K_m^*} - (\delta+\mu_m) I_m^v, \end{aligned}$$ \right. $

(14)

DFE of System (14) is $(S_f,S_m,I_f,I_f^v,I_m,I_m^v)=\bigg((1-p_f)K_f^*,(1-p_m)K_m^*,0,0,0,0\bigg)$. Using the next generation matrix method [17] we obtain the following two matrices

$ \mathcal{F}=\left[ $$\begin{smallmatrix} \sigma_fS_f\frac{I_m+I_m^v}{K_f^*+K_m^*} \\ (1-\tau) \sigma_f(K_f^*-S_f-I_f-I_f^v)\frac{I_m+I_m^v}{K_f^*+K_m^*}\\ \sigma_mS_m\frac{I_f+I_f^v}{K_f^*+K_m^*}\\ (1-\tau) \sigma_m(K_m^*-S_m-I_m-I_m^v)\frac{I_f+I_f^v}{K_f^*+K_m^*}\\ 0\\ 0 \end{smallmatrix}$$ \right]$

$\text{and}\;\; \mathcal{V}=\left[ $$\begin{smallmatrix} (\delta+\mu_f) I_f\\ (\delta+\mu_f) I_f^v\\ (\delta+\mu_m) I_m\\ (\delta+\mu_m) I_m^v\\ -r_f(1-p_f)K_f^*(1-\frac{K_f^*}{K_f}) + \sigma_fS_f\frac{I_m+I_m^v}{K_f^*+K_m^*} - \delta I_f + \mu_f S_f\\ -r_m(1-p_m)K_m^*(1-\frac{K_m^*}{K_m}) + \sigma_mS_m\frac{I_f+I_f^v}{K_f^*+K_m^*} - \delta I_m + \mu_m S_m \end{smallmatrix}$$ \right].$

Let $Q=(I_f,I_f^v,I_m,I_m^v ,S_f,S_m)=\bigg(0,0,0,0,(1-p_f)K_f^*,(1-p_m)K_m^*\bigg)$. Evaluating the Jacobian matrices of $\mathcal{F}$ and $\mathcal{V}$ at $Q$, we obtain the following matrices

$D\mathcal{F}(Q)=\left[ $$\begin{array}{c|c} F&0\\ \hline 0&0 \end{array}$$ \right]\;\;\text{and}\;\;D\mathcal{V}(Q)=\left[ $$\begin{array}{c|c} V&0\\ \hline W&U \end{array}$$ \right],$

where

$F=\left[ $$\begin{smallmatrix} 0&0&\frac{\sigma_f(1-p_f) K_f^*}{K_f^*+K_m^*}&\frac{\sigma_f(1-p_f) K_f^*}{K_f^*+K_m^*}\\ 0&0&\frac{(1-\tau)\sigma_fp_fK_f^*}{K_f^*+K_m^*}&\frac{(1-\tau)\sigma_fp_fK_f^*}{K_f^*+K_m^*}\\ \frac{\sigma_m(1-p_m) K_m^*}{K_f^*+K_m^*}&\frac{\sigma_m(1-p_m) K_m^*}{K_f^*+K_m^*}&0&0\\ \frac{(1-\tau)\sigma_mp_m K_m^*}{K_f^*+K_m^*}&\frac{(1-\tau)\sigma_mp_m K_m^*}{K_f^*+K_m^*}&0&0\\ \end{smallmatrix}$$ \right],$

$V=\left[ $$\begin{smallmatrix} \delta+\mu_f&0&0&0\\ 0&\delta+\mu_f&0&0\\ 0&0&\delta+\mu_m&0\\ 0&0&0&\delta+\mu_m \end{smallmatrix}$$ \right], U=\left[ $$\begin{smallmatrix} \mu_f&0\\ 0&\mu_m \end{smallmatrix}$$ \right]$

$ \text{and}\;\; W=\left[ $$\begin{smallmatrix} -\delta&0&\frac{\sigma_f(1-p_f) K_f^*}{K_f^*+K_m^*}&\frac{\sigma_f(1-p_f) K_f^*}{K_f^*+K_m^*}\\ \frac{\sigma_m(1-p_m) K_m^*}{K_f^*+K_m^*}&\frac{\sigma_m(1-p_m) K_m^*}{K_f^*+K_m^*}&-\delta&0 \end{smallmatrix}$$ \right].$

Hence,

$FV^{-1}= $$\begin{bmatrix} 0&0&\frac{\sigma_f(1-p_f)K_f^*}{(\delta+\mu_m)(K_f^*+K_m^*)}&\frac{\sigma_f(1-p_f)K_f^*}{(\delta+\mu_m)(K_f^*+K_m^*)}\\ 0&0&\frac{(1-\tau)\sigma_fp_fK_f^*}{(\delta+\mu_m)(K_f^*+K_m^*)}&\frac{(1-\tau)\sigma_fp_fK_f^*}{(\delta+\mu_m)(K_f^*+K_m^*)}\\ \frac{\sigma_m (1-p_m)K_m^*}{(\delta+\mu_f)(K_f^*+K_m^*)}&\frac{\sigma_m (1-p_m)K_m^*}{(\delta+\mu_f)(K_f^*+K_m^*)}&0&0 \\ \frac{(1-\tau)\sigma_m p_mK_m^*}{(\delta+\mu_f)(K_f^*+K_m^*)}&\frac{(1-\tau)\sigma_m p_mK_m^*}{(\delta+\mu_f)(K_f^*+K_m^*)}&0&0 \end{bmatrix}$$ .$

By the next generation matrix method [17], the reproduction number for Model (12), $\mathcal{R}_0^v$, is the spectral radius of the next generation matrix, $\rho(FV^{-1})$, and

$\mathcal{R}_0^v=\rho(FV^{-1})=\sqrt{\mathcal{R}_{0f}^v\mathcal{R}_{0m}^v},\text{where }$

$\mathcal{R}_{0f}^v=\frac{(1-\tau p_f)\sigma_f K_f^*}{(\delta+\mu_f)(K_f^*+K_m^*)}\;\;\text{and}\;\;\mathcal{R}_{0m}^v=\frac{(1-\tau p_m)\sigma_m K_m^*}{(\delta+\mu_m)(K_f^*+K_m^*)}.$

$\mathcal{R}_0^v<1$ implies DFE $((1-p_f)K_f^*,p_fK_f^*,(1-p_m)K_m^*,p_mK_m^*,0,0,0,0)$, is locally asymptotically stable, whereas $\mathcal{R}_0^v>1$ implies DFE is unstable [17].

$\mathcal{R}_{0f}^v=(1-\tau p_f)\mathcal{R}_{0f}\text{and}\mathcal{R}_{0m}^v=(1-\tau p_m)\mathcal{R}_{0m}.$

Since $0<1-\tau p_f<1$ and $0<1-\tau p_m<1$, $\mathcal{R}_{0}^v<\mathcal{R}_{0}<1$.

Thus, adopting a HPV vaccination program decreases the basic reproduction number, $\mathcal{R}_{0}$, in AA population.

Using the parameter values of Table 3, $\tau=90\%$ [14], $p_f=39\%$ [9] and $p_m=20.4\%$ [9], we obtain that $\mathcal{R}_{0f}^v=0.1243$ and $\mathcal{R}_{0m}^v=0.2346$. Hence, $\mathcal{R}_0^v=0.1708$.

In the next section, we use sensitivity analysis to illustrate the impact of model parameters on $\mathcal{R}_0^v$.

---

6.1. Sensitivity analysis of $\mathcal{R}_0^v$

We use Definition 4.1 to derive the sensitivity indices of the basic reproduction number $\mathcal{R}_0^v$ and we evaluate them using, $\tau=90\%$, $p_f=39\%$, $p_m=20.4\%$ and parameter values of Table (3).

From Table 7 and Figure 5, increasing (respectively, decreasing) the clearance rate, $\delta$, by $1\%$ will decrease (respectively, increase) the value of $\mathcal{R}_{0}^v$ by about $0.99\%$. Increasing (respectively, decreasing) the infection rate of the AAF population, $\sigma_f$, by $1\%$ will increase (respectively, decrease) the value of $\mathcal{R}_{0}^v$ by about $0.5\%$. Increasing (respectively, decreasing) the infection rate of the AAM population, $\sigma_m$, by $1\%$ will increase (respectively, decrease) the value of $\mathcal{R}_{0}^v$ by about $0.5\%$. Increasing (respectively, decreasing) the success rate of vaccination, $\tau$, by $1\%$ will decrease (respectively, increase) the value of $\mathcal{R}_{0}^v$ by about $0.38\%$. Increasing (respectively, decreasing) the proportion of HPV vaccinated females, $p_f$, by $1\%$ will decrease (respectively, increase) the value of $\mathcal{R}_{0}^v$ by about $0.27\%$.

---

6.2. African American male and female HPV model with vaccination simulations

To illustrate the impact of HPV on AAF and AAM populations of 16 years and older when a vaccination program is applied with $\tau=90\%$, $p_f=39\%$ and $p_m=20.4\%$, we simulate Model (12) with the parameter values listed in Table 3 and the initial conditions given in Table 6.

Simulations of our HPV Model (12) are performed using Matlab software, and are illustrated in Figure 6. Figures 6 (a-b) show that susceptible population of AAF of 16 years and older, $S_f$ (respectively, vaccinated AAF of 16 years and older $S_f^v$), increases over time. Similarly, Figures 6 (c-d) show that susceptible population of AAM of 16 years and older, $S_m$ (respectively, vaccinated AAM of 16 years and older $S_m^v$), increases over time. In the total population of African American population, $\mathcal{R}_{0}^v<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_f^v$), 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_m^v$), 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 $\tau=90\%$, while we vary $p_f$ and $p_m$.

Note that in Figure 6 and Table 6, $S_f(0)+S_f^v(0)$ (respectively, $I_f(0)+I_f^v(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_m^v(0)$ (respectively, $I_m(0)+I_m^v(0)$) is the same as $S_m(0)$ (respectivrly, $I_m(0)$) of Table 4 and Figure 4.

---

6.3. Impact of HPV vaccination

In AA population, $\mathcal{R}_0^v<\mathcal{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.

Furthermore, in both AAF and AAM populations, Figures 9 and 10 show that the number of infected populations is lower when the population is under a vaccination policy than when the population is not being vaccinated. Thus, HPV vaccines that provide partial immunity to both AAF and AAM populations of 16 years and older not only lower the number of HPV infectives but increase the number of susceptibles in both female and male populations.

Furthermore, we obtained in Figures 7-10 that the increase (respectively, decrease) in the susceptible (respectively, HPV infective) populations is larger when a bigger proportion of the population is vaccinated.

---

7. Conclusion

We use a two-sex HPV model with "fitted" logistic demographics to study HPV disease dynamics in AAF and AAM populations of 16 years and older. In agreement with Lee and Tameru [13], we obtained that in AA population, $\mathcal{R}_0<1$ and HPV cannot get started in a fully susceptible AA population.

Using sensitivity analysis on $\mathcal{R}_0$, we obtained the following results:

● Increasing (respectively, decreasing) the clearance rate, $\delta$, by $1\%$ will decrease (respectively, increase) the value of $\mathcal{R}_{0}$ by about $0.99\%$.

● Increasing (respectively, decreasing) the infection rate of the AAF population, $\sigma_f$, by $1\%$ will increase (respectively, decrease) the value of $\mathcal{R}_{0}$ by about $0.5\%$.

● Increasing (respectively, decreasing) the infection rate of the AAM population, $\sigma_m$, by $1\%$ will increase (respectively, decrease) the value of $\mathcal{R}_{0}$ by about $0.5\%$.

In the second part of the paper, we extended our model to include vaccination classes in both male and female AA populations of 16 years and older. We obtained that in AA population when the vaccination program is implemented, $\mathcal{R}_0^v<1$. Using sensitivity analysis on $\mathcal{R}_0^v$, we obtained the following results:

● Increasing (respectively, decreasing) the clearance rate, $\delta$, by $1\%$ will decrease (respectively, increase) the value of $\mathcal{R}_{0}$ by about $0.99\%$.

● Increasing (respectively, decreasing) the infection rate of the AAF population, $\sigma_f$, by $1\%$ will increase (respectively, decrease) the value of $\mathcal{R}_{0}$ by about $0.5\%$.

● Increasing (respectively, decreasing) the infection rate of the AAM population, $\sigma_m$, by $1\%$ will increase (respectively, decrease) the value of $\mathcal{R}_{0}$ by about $0.5\%$.

● Increasing (respectively, decreasing) the success rate of vaccination, $\tau$, by $1\%$ will decrease (respectively, increase) the value of $\mathcal{R}_{0}$ by about $0.38\%$.

● Increasing (respectively, decreasing) the proportion of HPV vaccinated females, $p_f$, by $1\%$ will decrease (respectively, increase) the value of $\mathcal{R}_{0}$ by about $0.27\%$.

Furthermore, using the extended model with vaccination we obtained the following results:

● Adopting a vaccination policy lowers HPV infections in both AAF and AAM populations.

● Vaccinating a larger proportion of AAF and AAM populations leads to fewer cases of HPV infections in the vaccinated population.

---

Acknowledgments

This research was partially supported by National Science Foundation under grant DUE-1439758.

---