Comparative analysis of practical identifiability methods for an SEIR model

1.   Introduction

Mathematical models are used to study sexually transmitted infections (STIs) and to inform policy, practice, and resource allocation. For example, the models can be used to study how STIs spread, the impact of interventions, and the potential effect of vaccines. When modeling STIs, some important factors to consider are: Duration and type of partnerships (for example, monogamous partnerships can protect susceptible individuals from infection) ^[1]; granularity of the model (enabling capturing patterns of sexual contact at different levels, such as individual, partnership, or network level) ^[2]; structuring of the population (in groups vs. explicitly simulating individuals and their contacts, which provides more detail and can include microstructures in the sexual partnership network) ^[3]. The role of casual partners, which can be enabled by online dating applications, is also becoming increasingly important ^[4].

Research has shown that sexual contact patterns differ among individuals according to sexual orientation, both in terms in actual behaviors and in frequency and number of sexual partners. For example, Blanc et al. (2023) ^[5] conducted a study of 2288 young adults that reported

"Bisexual women engaged in a greater number of sexual behaviors and had more positive attitudes toward sexual behaviors than heterosexual women. Homosexual men engaged in a greater number of sexual behaviors than heterosexual men, and homosexual and bisexual men had more positive attitudes toward sexual behaviors than heterosexuals. Finally, we show the mediating role of attitudes toward sexual behaviors in the relationship between sexual orientation and the number of sexual behaviors engaged in."

Therefore, to understand the dynamics of transmission of sexually transmitted infections, it is important to consider these differences that lead to non-homogeneous transmission rates across sexes and sexual orientations. It may be desirable to divide the sexually active population by sex (female vs. male), and within each of the sexes further subdivide by sexual orientation into heterosexual (individuals who have sexual contact exclusively with those of the opposite sex), homosexual (individuals who have sexual contact exclusively with those of the same sex), and bisexual (individuals who have sexual contact with those of both sexes). However, such structure would require the identification of 8 effective contact rates for each of the sexes: Four affecting bisexual individuals –those with homosexual and bisexual individuals of the same sex, and those with heterosexual and bisexual individuals of the opposite sex; two affecting heterosexual individuals –those with heterosexual and bisexual individuals of the opposite sex; and two affecting homosexual individuals –those with homosexual and bisexual individuals of the same sex.

For the practical use of such a model, one needs to identify far too many parameters –including the sixteen effective per capita transmission rates– from too few data points, because for most sexually transmitted infections there is almost no data on prevalence and very scarce data on incidence.

To circumvent this serious barrier, we take the approach of ^[4] and consider only two subpopulations instead of six, breaking up the total sexually active population into two sexes, female and male, and summarizing all sexual contacts by including average per capita effective transmission rates across the only four possibles routes across the sexes. This is a limitation in the model because it is equivalent to assuming every individual in the population is bisexual and has the prescribed mean transmission rate towards each of the sexes. However, this approach significantly reduces the parameter identification problem, making it manageable, so that we may estimate the model parameters from data more reliably.

In contrast with the SIS model proposed and analyzed in ^[4], which was addressing bacterial infections with treatments that do not confer immunity, we propose and analyze here a SIRS model for viral infections with treatments that do confer immunity, be it permanent or temporary. We have also made a new estimation of the sexually active population by sex and age that significantly improves that of ^[6].

2.   Mathematical model

We present a mathematical framework to capture the dynamics of viral, sexually transmitted infections, such as HIV, including partner selection across two sexes. For simplicity, we do not consider age-structured populations but rather structure the sexually active population only by sex. The transmission of HIV, for example, is mainly through sexual interaction ^[7]. We consider two sex groups, female and male, and a SIRS epidemiological model with sexual interactions among them in all combinations. With minor modifications, the model can be applied to any viral STI. The total population for each group is the sexually active population of the corresponding sex, with sizes denoted as $N_f = N_f(t)$ and $N_m = N_m(t)$, respectively. Because we are interested in long-term dynamics, the population in each compartment will be subject to their natural recruitment, mortality, and emigration. We subdivide the population of females ($k = f$) and males $(k = m)$ into six mutually exclusive compartments, three for each sex: Susceptible (of size $S_k$), infected (of size $I_k$), and recovered individuals (of size $R_k$).

Susceptible ($S_k$)

The populations of susceptible females and males grow from the recruitment of newly sexually active individuals at a rate $\Lambda_k$, and from previously immune individuals who lose their immunity to HIV-infection. They decrease at a per capita rate $\lambda_g$ because of infection, at a per capita rate $\mu_g$ because of cessation of sexual activity for any reason—including natural death.

To make the model more tractable in terms of the number of state variables (subpopulations) and parameters it requires, we introduce a limitation to our model by combining all sexual interactions with individuals of the same sex and of the opposite sex into a single force of infection, which is equivalent to assuming that all individuals are bisexual. Therefore, the force of infection for the female population is given by

where $\beta_{mf}$ is the per capita transmission rate from infected male to uninfected female and $\beta_{ff}$ is the per capita transmission rate from infected female to uninfected female.

Similarly, for the male population, the force of infection is given by

where $\beta_{mm}$ is the per capita transmission rate from infected male to uninfected male, and $\beta_{fm}$ is the per capita transmission rate from infected female to uninfected male. We assume constant per capita contact rates for all combinations of sexes.

Infected ($I_k$)

For simplicity, we do not consider a latency or incubation period, and infected individuals are identified with the infectious in our model, whether or not they present clinical symptoms or are detected/confirmed by HIV test, Pap test, PCR tests, etc. For each sex, the number of infected increases because of infection and decreases by clearance of the virus—at per capita rate $\gamma_{k}$, by cessation of sexual activity for any reason other than infection-induced death—at per capita rate $\mu_k$, and by infection-induced death—at per capita rate $\delta_k$.

Recovered/Immune ($R_k$)

For each sex, this population grows due to the clearance of the infection with temporary immunity, and decreases from cessation of sexual activity for any reason—at per capita rate $\mu_k$, and from loss of immunity at per capita rate $\theta_k$.

Two-sex SIRS model

Our model is given by

where the total size of the population of sex $k$ is given by

solution of the initial value problem

It follows that $N_k$ is bounded above by its obvious super-solution $\bar N_k$ corresponding to $\delta_k = 0$,

We show in Figure 1 a flow diagram of the model, where all the parameters included therein are per capita rates, except for the recruitment rates $\Lambda_k$.

The summary of all parameters in the model and their description are in Table 1.

3.   Mathematical analysis

Lemma 3.1. The closed set

is positively-invariant and attracting for the flow of model (2.3). Here,

are the asymptotic values of the supersolutions (2.6), $\bar N_k$, $k = f, m$.

The lemma proves that the model is biologically meaningful in that it preserves non-negativity of solutions and their boundedness by the most natural upper bound, which is the total population size in the absence of infection for each sex.

3.1.   Disease-free equilibrium and basic reproduction number

Equating the sizes of the infected classes $I_f$ and $I_m$ and the right-hand sides of all equations in (2.3) to zero, we see that the model has the disease-free equilibrium (DFE) given by

where

Definition 3.1 (Basic Reproduction Number). The basic reproduction number, $\mathscr{R}_0$, for the model of Figure 1 is the expected number of secondary infections generated by a single infected individual during her/his entire infectious period when introduced into a completely susceptible population.

The basic reproduction number for the female-only SIRS model is

and the basic reproduction number for the male-only SIRS model is

The numerators in (3.4) and (3.5) are the effective per capita transmission rates and the denominators are the per capita loss of infectivity rates for each sex. Thus, the fraction equals the product of the per capita production rate of new infections times the expected value of the infectious period, which exactly matches the definition.

For the two-sex model (2.3), it can be computed as the spectral radius of the next generation matrix (details in the Appendix)

where

and

Because $\mathscr{R}_0$ is the spectral radius of the positive operator $\mathcal{G}$, it is equal to its largest eigenvalue,

where,

and

We can prove some useful relations between the three basic reproduction numbers (3.4)–(3.6) (see Appendix).

Proposition 3.1.

The first statement in part a. of Proposition 3.1 expresses the intuitive idea that the two-sex transmission model should have a larger basic reproduction number than either one-sex transmission model because the intersex transmission can only increase the number of secondary cases from a primary case under just intrasex transmissions. The second statement just says that such an increase will not exist only when intersex transmission is not present.

The first statement of part b. expresses the idea that when the dominant transmission pathways are intersex, the introduction of a single infected individual into the two-sex population will produce at least as many secondary cases as the introduction of one primary case of each sex resulting in the total number of secondary cases exceeding the sum of the secondary cases from the two intrasex transmission pathways. The second statement just says that such excess will not exist when intersex transmissions are balanced with intrasex transmissions, in which case the total number of secondary cases will equal the sum of the secondary cases from the two intrasex transmission pathways.

For the case of identical one-sex basic reproduction numbers, we have the following:

Corollary 3.1.1. If $\mathscr{R}_{0, f} = \mathscr{R}_{0, m} = {\mathscr{R}_{0g}}$, then

This corollary quantifies how much larger the basic reproduction number (3.6) is than the one-sex basic reproduction numbers when the latter are equal. It is a special case of part a. of Proposition 3.1 that makes the increase precise in terms of the ratio of intersex to intrasex transmission rates.

We now turn to the local asymptotic stability of the disease-free equilibrium. Its local asymptotic stability is controlled by the basic reproduction number.

Proposition 3.2. [Stability of DFE] $E_0^*$ is locally-asymptotically stable (LAS) if ${\mathscr{R}_0} < 1$ and unstable if ${\mathscr{R}_0} > 1.$

The proof (detailed in the Appendix) is based on the linearization of the model around $E_0^*$ and proving that a sufficient condition for all the eigenvalues of the Jacobian matrix at $E_0^*$ to have negative real part is precisely ${\mathscr{R}_0} < 1$.

3.2.   Endemic equilibria

Model (2.3) has two nontrivial boundary ("only-one-sex-infected") equilibria, i.e., those where $I_f^* = 0$ and $I_m^* > 0$, or $I_m^* = 0$ and $I_f^* > 0$. The "Only-Males-Infected" Endemic Equilibrium, with $I_m^* > 0$ and $I_f^* = 0$ (equivalent to $\beta_{mf} = 0$ and $\beta_{mm} > 0$) is ${E^*_m = \left(N^*_f, 0, 0, S^*_m, I^*_m, R^*_m\right)}$, given by

The "Only-Females-Infected" Endemic Equilibrium, with $I_m^* = 0$ and $I_f^* > 0$ (equivalent to $\beta_{fm} = 0$ and $\beta_{ff} > 0$), is ${E^*_f = \left(S^*_f, I^*_f, R^*_f, N^*_m, 0, 0\right)\!}$ given by

It is very easy to verify that (3.7) and (3.8) indeed satisfy $S^*_m+I^*_m+R^*_m = N_m^*$ and $S^*_f+I^*_f+R^*_f = N_f^*$, and we have the following results regarding their existence and local asymptotic stability.

Proposition 3.3. The "only-males-infected" equilibrium $E^*_m$ given by (3.7) exists in $\mathscr{D}$ and differs from the DFE if, and only if, ${\mathscr{R}}_{0, m} > 1$. It is locally asymptotically stable if, and only if, ${\mathscr{R}}_{0, m} > 1$ and ${\mathscr{R}}_{0, f} < 1$.

Proposition 3.4. The "only-females-infected" equilibrium $E^*_f$ given by (3.8) exists in $\mathscr{D}$ and differs from the DFE if, and only if, ${\mathscr{R}}_{0, f} > 1$. It is locally asymptotically stable if, and only if, ${\mathscr{R}}_{0, f} > 1$ and ${\mathscr{R}}_{0, m} < 1$.

Identification of all equilibria in the interior of the conic domain $\mathscr{D}$ is a bit more complicated. In the case of only one intersex transmission pathway (i.e., $\beta_{fm} = 0$ or $\beta_{mf} = 0$), the unique interior equilibrium for each case can be identified explicitly. We have

Proposition 3.5. Let $\mathscr{R}_{0, f} > 1$ and $\mathscr{R}_{0, m} > 1$ —whereby also $\mathscr{R}_0 > 1$ by Proposition 3.1. If $\beta_{fm} = 0$ (i.e., females cannot infect males), then the model admits a unique endemic equilibrium ${E^* = (S^*_f, I^*_f, R^*_f, S^*_m, I^*_m, R^*_m)\!}$, given by

Proposition 3.6. Let $\mathscr{R}_{0, f} > 1$ and $\mathscr{R}_{0, m} > 1$ —whereby $\mathscr{R}_0 > 1$ by Proposition 3.1. If $\beta_{mf} = 0$ (i.e., males cannot infect females), then the model admits a unique endemic equilibrium ${E^* = (S^*_f, I^*_f, R^*_f, S^*_m, I^*_m, R^*_m)\!}$, given by

The quantities $a_i$ and $b_i$, $1\le i\le4$ are given by (5.29) and (5.32), and are the coefficients of the quadratic system (5.33) for equilibria of $I_m$ and $I_f$ (see the Appendix).

For the general case when all transmission pathways are available, we do not have explicit expressions for interior equilibria but can describe them completely in terms of up to three positive $I_m^*$ values given as roots of the cubic polynomial $\alpha_0+\alpha_1I_m+\alpha_2I_m^2+\alpha_3I_m^3$ with coefficients (5.38)–(5.41) that contain derived parameters $\sigma_0$, $\sigma_1$, and $\sigma_2$ given by (5.42).

By Descartes' Rule of Signs (see ^[8]), the maximum number of positive real roots a polynomial may have is given by the number of sign changes in consecutive pairs of its coefficients. Because (5.41) implies $\alpha_3 < 0$, the only way for 3 distinct positive roots to exist is that $\alpha_0 > 0$, $\alpha_1 < 0$, and $\alpha_2 > 0$. Moreover, the only way for no positive real roots to exist is that $\alpha_0, \alpha_1, \alpha_2 < 0$, that is, $\sigma_0, \sigma_1, \sigma_2 < 1$.

We summarize here the results we obtain for interior equilibria.

Proposition 3.7. [Existence and stability of of endemic equilibria] Endemic equilibria for model (2.3) exist in the following cases:

1) The boundary equilibrium $E^*_{m} = (S_m^*, I_m^*, R_m^*, S_f^*, 0, 0)$ exists and is locally asymptotically stable if $\mathscr{R}_{0, m} > 1$ and $\mathscr{R}_{0, f} < 1$.

2) The boundary equilibrium $E^*_{f} = (S_m^*, 0, 0, S_f^*, I_f^*, R_f^*)$ exists and is locally asymptotically stable if $\mathscr{R}_{0, m} < 1$ and $\mathscr{R}_{0, f} > 1$.

3) At least one positive equilibrium $E^* = (S_m^*, I_m^*, R_m^*, S_f^*, I_f^*, R_f^*)$ exists if at least one of $\sigma_0, \sigma_1$, or $\sigma_2$ is bigger than 1, and at least one of the positive roots $I^*_m$ of (5.43) satisfies the bounds given in (5.37).

Remark. We see in Table 3 in the Appendix that a necessary condition for a unique positive equilibrium to exist is

Remark. The classical threshold condition for the existence of a positive equilibrium (at least one) is $\mathscr{R}_0 > 1$, and we can see that when there is no disease-induced mortality ($\delta_f = \delta_m = 0$),

even though $\mathscr{R}_{0}\ge\max\{\mathscr{R}_{0, f}, \mathscr{R}_{0, m}\} > 1$. In fact, (3.6) and (5.16) imply that

whereby

may even vanish. Sufficient conditions for $\sigma_0$ to be greater than 1 are $\mathscr{R}_{0, f} > 1, \mathscr{R}_{0, m} > 1$, and $\beta_{ff}\, \beta_{mm}\le\beta_{fm}\, \beta_{mf}$ because then, by (3.1),

4.   Numerical simulations

We consider the sexually active population of the U.S. as described in ^[6], consisting of individuals aged at least 15 who have engaged in sexual contact (oral, vaginal, or anal) with other individuals of any sex.

The annual numbers of sexually active individuals by sex in the United States is estimated using the annual sex-and-age structure of the U.S. population ^[9] and multiplying each cohort's size by the estimated percentage of sexually active within the corresponding sex-and-age cohort. Those percentages are taken from ^[10,11] for ages 15–20 (using data from 2011 to 2015 and from 2015 to 2017), from ^[12] for ages 21–29, from ^[13] for ages 30–39 (using data for 2006–2008), and from ^[14] for ages 40–44 (using data for 2006–2008). For people 50 years old and older, we use the data from the meta study about sexual activity among older adults ^[15]. The sizes of all sexually active age cohorts for each sex are then summed to produce the estimated number of sexually active men and women.

The percentages of newly sexually active individuals by age and sex are computed as the difference of the percentages of sexually active individuals for each age group and sex, and the percentage of sexually active individuals for that sex in the preceding age group. Their sums over all age cohorts for each sex provide an estimate for the rate of newly sexually active men and women entering the sexually active population ($\Lambda_m$ and $\Lambda_f$, respectively).

For the purpose of forecasting, we subsequently model the sizes of the sexually active populations of each sex through the following "learning curve" functions:

where the initial sizes of the sexually populations of men and women, and four parameters therein ($\Lambda_g\! =$ rate at which sexually inactive people of sex $g$ enter the sexually active population, $\mu_g =$ per capita rate at which sexually active individuals of sex $g$ exit the sexually active population, $g = f, m$) have the following values, obtained by fitting the model (4.1) and (4.2) to the estimated annual numbers of sexually active men and women as described above (population numbers and rates are in 1000's):

We depict in Figure 2 the resulting graphs of (4.1) and (4.2) together with the estimated sexually active population by sex for the years 2009–2019, showing an excellent fit of model to data.

We next show the results of simulations corresponding to different parameter values that result in various combinations of signs for $\mathscr{R}_0-1$, $\mathscr{R}_{0, f}-1$, and $\mathscr{R}_{0, m}-1$, with consequent different asymptotic behaviors represented in Propositions 3.2–3.7. For all our simulations, we assum no infection-induced mortality ($\delta_f = \delta_m = 0$) and use the demographic parameters from (4.3), fixing the per capita recovery rates and loss of immunity rates as

and initialized the infection SIRS model with the following cohort sizes:

Case 1: $\mathscr{R}_0 < 1$. According to Proposition 3.1, this entails $\mathscr{R}_{0, f} < 1$ and $\mathscr{R}_{0, m} < 1$. We set the effective per capita infectious contact rates as follows:

These result in $\mathscr{R}_0 = 0.98$, $\mathscr{R}_{0, f} = 0.15$, and $\mathscr{R}_{0, m} = 0.65$, which correspond to the DFE being locally asymptotically stable according to proposition 3.2. Because it is the only equilibrium point of the system, it turns out that it is globally asymptotically stable. We see in Figure 3 that the infected cohorts of both sexes disappear (left panel), as do those of recovered individuals (right panel). We also see in the right panel that the susceptible populations of men and women approach their limiting values $N_m^* = 99, \!267, \!000$ and $N^*_f = 89, \!433, \!000$, respectively.

Case 2: We changed only the effective per capita infectious contact rate from men to men slightly, to $\beta_{mm} = 0.53$ while keeping all other parameters the same as in Case 1:

These result in $\mathscr{R}_0 = 1.002 > 1$, while $\mathscr{R}_{0, f} = 0.15 < 1$, and $\mathscr{R}_{0, m} = 0.68 < 1$, destabilizing the DFE according to Proposition 3.2 while a unique endemic equilibrium bifurcates as $\mathscr{R}_{0}$ crosses the value 1 from left to right (according to Proposition 3.7). In this case, there are no boundary equilibria ("only-one-sex-infected") according to Proposition 3.7. The positive equilibrium –for which we do not have an explicit formula– is locally asymptotically stable, as shown in Figure 4.

Note that the results of Propositions 3.3 and 3.4 are equivalent, the first resulting in the existence and local asymptotical stability of $E^*_m$ and the second of $E^*_f$. We present simulations to exemplify just one of them.

Case 3: We set the effective per capita infectious contact rates as follows:

These result in $\mathscr{R}_0 = 1.5$, $\mathscr{R}_{0, f} = 1.5$, and $\mathscr{R}_{0, m} = 0.45$, corresponding to the "only-females-infected" boundary equilibrium being asymptotically stable according to Proposition 3.4, as shown in Figure 5.

Note that the results of Propositions 3.5 and 3.6 are analogous to each other and they correspond to endemic equilibria that are LAS and for which we have explicit formulas, the first having no transmission from woman to man and the second no transmission from man to woman. We present simulations to exemplify just one of them.

Case 4: We set the effective per capita infectious contact rates as follows:

These result in $\mathscr{R}_0 = 1.25$, $\mathscr{R}_{0, f} = 1.25$, and $\mathscr{R}_{0, m} = 1.1$, corresponding to the endemic equilibrium being asymptotically stable according to Proposition 3.6, as shown in Figure 6.

4.1.   Estimation of per capita transmission rates

Because there is neither recovery from nor immunity to HIV infection, we need to ensure that the compartments of recovered individuals are empty, that is

We accomplish this by selecting vanishing per capita recovery and loss-of-immunity rates:

Then, we use model (2.3) to fit the per capita transmissibility parameters $\beta$ to incidence data for HIV for the years 2008–2013 and obtained the following values:

The fitting is done by minimizing the maximum of the relative $l^\infty$-errors in incidence among men and women.

As we might expect, trasmissibility from woman to woman is two orders of magnitude smaller than those for any other sex combinations. We also see that transmissibility from man to woman seems to be about 20% smaller than from man to man, and that transmissibility from woman to man seems to be the largest.

We can see in Figure 7 and in Table 2 a very good fit of simulated incidence —represented by $\lambda S$ for each of the sexes— to actual incidence data, with errors for every year remaining below ±0.7% for males and below ±5.63% for women. The disparity in relative errors across sexes is due to the four-to-fivefold incidence among men as comapred to women.

4.2.   Model-based prediction

Using the best-fit parameters from (4.10), we use the model to forecast male and female incidence for the following two years. We see in the second panel of Figure 8 that the model-predicted values for incidence in women closely match the actual data, with errors just at 0.3% for 2014 and 1% for 2015.

The excellent predictive value of the model in this case is easily explained because the trend line we use for the years 2008–2013 for female incidence (the quadratic regression curve of the 6 annual data points) remain very consistent with the actual incidence in the following two years, as seen from the data points for 2014 and 2015 being very close to the quadratic regression curve.

In the first panel of Figure 8, we see a very different predictive quality: The forecast from the model values significantly underestimate the real data, with relative errors of about 5.2% for 2014 and 7% for 2015. This decreased predictive value is also easily explained by the shift in the trend of HIV incidence for men after 2013, which the model cannot account for because it assumes consistent transmission dynamics over time, represented by the cubic regression curve of the incidence data for men for the period from 2008 to 2013, represented by the dotted line.

4.3.   Sensitivity analysis

To better understand the impact that each of the four types of sexual contacts have on the basic reproduction number, we perform sensitivity analysis using the normalized sensitivity indexes of $\mathscr{R}_0$ with respect to $\beta_{ij}$, for $i,j\in\{f,m\}$.

Definition 4.1. [Normalized Sensitivity Index] The normalized sensitivity index, $SI_\rho(U)$, of a quantity of interest $U$ with respect to the parameter $\rho$, for $U\not =0\text{ and }\rho\not=0,$ is defined as

In Figure 9, we depict the normalized sensitivity indexes of $\mathscr{R}_0$ with respect to all behavioral parameters $\beta_{i,j}$ with $i,j\in\{m,f\}$, given by (4.11). To construct this graph, we set $\mu_m=0.034$ and $\mu_f=0.044$ from (4.3), and $\delta_m=\delta_f=0.05$ corresponding to the estimated HIV-induced per capita mortality rate, ^[16]. The per capita effective transmission rates $\beta_{i,j}$ are set at their baseline values identified in the parameter fitting described above, with the exception of the one for which we are computing the sensitivity.

It is interesting to see in Figure 9 that, even though the values of the effective per capita transmission rates from woman to man and from man to man are very similar, the basic reproduction number is over 4 times more sensitive to the latter than to the former. Also interesting is the fact that the sensitivity index of $\mathscr{R}_0$ with respect to $\beta_{ff}$ is increasing at an increasing rate, while those with respect to the other three effective per capita transmission rates are increasing at decreasing rates.

Similarly, we find the normalized sensitivity indices of $I_m^*$ and $I_f^*$ with respect to the three nonzero effective per capita transmission rates, for the endemic equilibrium values of Propositions 3.5 and 3.6, respectively. The changes made in parameter values are $\delta_m=\delta_f=0$, required to have the explicit expressions for $I_m^*$ and $I_f^*$ in these Propositions, and baseline values for the behavioral parameters $\beta_{i,j}$ with $i,j\in\{m,f\}$ such that Propositions 3.5 and 3.6 hold. Note that the necessary condition for the existence of the endemic equilibrium in Proposition 3.5, $\mathscr{R}_{0,f}>1$, is equivalent to $\beta_{ff}>\gamma_f+\mu_f=0.694$, while that for the existence of the endemic equilibrium in Proposition 3.6, $\mathscr{R}_{0,m}>1$, is equivalent to $\beta_{mm}>\gamma_f+\mu_f=0.784$. The corresponding sensitivity curves are graphed in Figure 10, where the curve corresponding to $\beta_{ff}$ on the left panel (in red) should be ignored for values of $\beta_{ff}\le0.694$, that is to the left of the vertical asymptote displayed, and the curve corresponding to $\beta_{mm}$ on the right panel (in black) should be ignored for values of $\beta_{ff}\le0.784$, that is to the left of the vertical asymptote displayed.

It is interesting to note in both panels that, for moderately large values of the effective per capita infection-transmission rates, the sensitivity of the equilibrium number of infected is pretty much the same with respect to any of the three nonzero $\beta$s, and quite small. Also noteworthy, is the fact that for the sex that can infect the opposite sex, the effective per capita same-sex transmission rate has a critical value at which the relative sensitivity index of its equilibrium value of infected switches from negative to positive. This means that for the effective per capita same-sex transmission rate that is small enough, a slight increase in its value is beneficial in reducing the size of the endemic equilibrium of infected within that sex and that size reaches a minimum at the critical value of the effective per capita same-sex transmission rate.

5.   Conclusions

We have described and analyzed a model for transmission of sexually transmitted infections across the four sexual routes (man-to-woman, woman-to-man, man-to-man, and woman-to-woman), not separating individuals by sexual preference (heterosexual, homosexual, and bisexual) because of the large increase in the number of parameters such a model requires. We kept all sexual contacts for the female and male populations occurring both with women and men, resulting in the undesirable implicit assumption that everyone in the population is bisexual. This results, for example, in needing only 4 effective per capita transmission rates (from infected men and women to uninfected men and uninfected women) rather than 16 that are needed for all possible routes of transmission in the sexual preference model (for each of the sexes, two for heterosexual compartments, two for homosexual compartments, and four for bisexual ones).

We have shown how to model the sexually active population of each sex using a learning curve model and, based on this model, how we may describe the transmission dynamics of viral STIs (such as HIV, herpes simplex type 2, and HPV) through a two-sex SIRS ODE model.

If we built a more complex demographic model, structured by sexual preferences as described above, possibly also separating ages 15–24 where almost half of the incidence occurs, we expect to find an improved fitting of incidence data and an added benefit of being able to compare effective transmission rates among sexual preference groups.

Our model may have up to 6 equilibrium points, the infection-free and two other boundary points of the positive orthant —corresponding to only-one-sex-infected, and up to 3 endemic equilibrium points in the interior of the positive orthant. We defined three basic reproduction numbers of the infection —one for each sex, $\mathscr{R}_{0,f}$ and $\mathscr{R}_{0,m}$, and one for the two-sex interactions, $\mathscr{R}_{0}$— that lead to the classical threshold conditions for long term dynamics depending on which are smaller or greater than 1. We have shown that, when $\mathscr{R}_0 <1$ (and a fortiori $\mathscr{R}_{0,f}<1$ and $\mathscr{R}_{0,m}<1$), the only equilibrium is the disease-free and it is locally asymptotically stable, and when $\mathscr{R}_{0,f}>1$ the only-females-infected equilibrium appears and is locally asymptotically stable if $\mathscr{R}_{0,m}<1$ and unstable if $\mathscr{R}_{0,m}>1$. Analogously, when $\mathscr{R}_{0,m}>1$ the only-males-infected equilibrium appears and is locally asymptotically stable if $\mathscr{R}_{0,f}<1$ and unstable if $\mathscr{R}_{0,f}>1$. The disease-free equilibrium always exists, while the only-one-sex-infected boundary equilibrium points bifurcate from it when the corresponding single-sex basic reproduction numbers cross the value 1 from left to right. When both do, both boundary equilibria lose their local asymptotic stability, and up to three endemic equilibria may appear.

The model provides a natural way to estimate prevalence of the STI under consideration, which is very useful because for most STIs, prevalence data is very scarce, while incidence data is usually more readily available. Using just one point in time as the initial time for which prevalence is known, the parameter fitting can be done by minimizing the error function chosen for differences between incidence data and estimated incidence from the model given by $\lambda S$.

A more detailed study of endemic equilibrium points and their stability would be useful, particularly to understand if it is possible to have more than one with local asymptotic stability and what their basins of attraction are.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare there is no conflict of interest.

Appendix

We provide details of the proofs of some theoretical results presented in the Section "Mathematical Analysis".

Lemma 3.1.

The closed set

is positively-invariant and attracting for the flow of model (2.3). Here,

are the asymptotic values of the supersolutions (2.6), $\bar N_k$, $k = f, m$.

Proof. First, let us recall that, for $k = f, m$,

Thus, using (5.2), we have

Therefore, all forward trajectories for our model (2.3) are bounded.

Now, let us prove that for non-negative initial conditions the trajectories stay in the first orthant.

If $I_f(0) = I_m(0) = 0$, then $I_k,$ and $R_k$ vanish for all time, and $S_k$ will monotonically approach $N_k^*$. To prevent this trivial dynamics, we assume that $I_f(0)\not = 0$ or $I_m(0)\not = 0$. Then, $\lambda_f(t) > 0$ and $\lambda_m(t) > 0$ for $t > 0$, which is sufficiently small. Hence, all six state variables of (2.3) are positive for $t > 0$ sufficiently small. Now, let

We want to prove that $t^\prime = +\infty$. Suppose it is finite. The first and fourth equations of model (2.3) are ($g = f, m$)

whereby

Similarly, for $k, l\in\{f, m\}$,

so that

Finally, following the same reasoning for the last equation,

leading to

Hence, by (5.4)–(5.7), we have proved that $\bar t = +\infty$ and all trajectories beginning in the positive orthant stay in its interior for all time. It now follows from (5.2) that $\mathscr{D}$ defined in (5.1) is positively invariant.

To compute the basic reproduction number, we proceed as in ^[17], using the next generation matrix. The matrix $\mathcal{F}$ containing the unit rates of new infections and the matrix $\mathcal{V}$ containing the unit rates of transitions out of infected classes associated with our model are given by

where

and

where,

The next generation matrix, $\mathcal{G}$, is

Lemma 5.1. The basic reproduction number is given by

where,

with

the basic reproduction number for the female-only SIRS model,

the basic reproduction number for the male-only SIRS model, and

Proof. Note that, by definition, we have $\mathscr{R}_0 = \rho(\mathcal{F} \mathcal{V}^{-1})$, where $\rho(\; )$ represents the spectral radius of the positive operator $\mathcal{G}$, which is its largest eigenvalue. Moreover,

Since all the parameters are non-negative, it follows from (5.12) that $\tau > 0$. We can verify that the radicand in (5.11) is non-negative, thus verifying $\mathscr{R}_0$ is a positive real. We can also see that the radicand vanishes if, and only if, there is no two-way heterosexual transmission ($\beta_{fm} = 0$ or $\beta_{mf} = 0$) and $\mathscr{R}_{0, f} = \mathscr{R}_{0, m}$.

Indeed,

We now prove the following relations between the three basic reproduction numbers.

Proposition 3.1. For the three basic reproduction numbers (5.11), (5.13) and (5.14), the following relations hold:

a. $\mathscr{R}_{0}\ge\max\left\{\mathscr{R}_{0, f}, \mathscr{R}_{0, m}\right\}, \quad \mbox{and}\quad \mathscr{R}_{0} = \max\left\{\mathscr{R}_{0, f}, \mathscr{R}_{0, m}\right\}\iff\beta_{fm}\beta_{mf} = 0$.

That is, the basic reproduction number $\mathscr{R}_{0}$ is at least as large as the largest of the single-sex reproduction numbers $\mathscr{R}_{0, f}$ and $\mathscr{R}_{0, m}$ with equality holding if, and only if, the product of the opposite-sex per capita transmission rates vanishes (that is, there is no two-way heterosexual transmission);

b. $\beta_{ff}\beta_{mm}\le\beta_{fm}\beta_{mf}\Longrightarrow \mathscr{R}_{0}\ge\mathscr{R}_{0, f}+\mathscr{R}_{0, m}, \quad \mbox{and}\quad \mathscr{R}_{0} = \mathscr{R}_{0, f}+\mathscr{R}_{0, m}\iff\beta_{ff}\beta_{mm} = \beta_{fm}\beta_{mf}$.

That is, if the product of the same-sex per capita transmission rates is no larger than the product of the opposite-sex ones, then $\mathscr{R}_{0}$ is at least as large as the sum of the single-sex reproduction numbers $\mathscr{R}_{0, f}+\mathscr{R}_{0, m}$ with equality holding if, and only if, the product of the same-sex per capita transmission rates is equal to the product of the opposite-sex ones.

Proof. Note that (5.16) can be rewritten as $\tau^2-4\mathcal{D} = \left(\mathscr{R}_{0, f}-\mathscr{R}_{0, m}\right)^2+\frac{4\beta_{fm}\beta_{mf}}{\Upsilon_f\Upsilon_m}$ which implies that

with equality holding if, and only if, one of the opposite-sex per capita transmission rates vanishes.

The relations in part a. of the proposition now follow from (5.11) and (5.17):

with equality holding if, and only if, $\beta_{fm}\beta_{mf} = 0$, as needed.

The relations in part b. of the proposition follow from (5.11) and (5.15) because $\beta_{ff}\beta_{mm}\le\beta_{fm}\beta_{mf}$ is equivalent to $\mathcal{D}\le0$ and $\beta_{ff}\beta_{mm} = \beta_{fm}\beta_{mf}$ to $\mathcal{D} = 0$.

We turn to the case of identical one-sex basic reproduction numbers.

Corollary 3.1.1. If $\mathscr{R}_{0, f} = \mathscr{R}_{0, m} = {\mathscr{R}_{0g}}$, then

Proof. The relation follows from (5.11) and (5.16).

We now prove a classic result on the local asymptotical stability of the disease-free equilibrium.

Proposition 3.2. $E_0^*$ is locally-asymptotically stable (LAS) if ${\mathscr{R}_0} < 1$ and unstable if ${\mathscr{R}_0} > 1.$

Proof. We follow the argument of ^[4]. The Jacobian of the vector function of right-hand sides of our model's flow at the DFE is

where $\omega_{fm}$ and $\omega_{mf}$ are given by (5.8), and $\Upsilon_f$, $\Upsilon_m$ are given by (5.9), $\mathscr{R}_{0, f}$ is given by (5.13), and $\mathscr{R}_{0, m}$ is given by (5.14).

The eigenvalues of $J(E^*_0)$ are

and the eigenvalues of the matrix

namely

We see that these are both real, with $\lambda_5 < \lambda_6$. They are both negative if, and only if

The second condition in (5.19) is

equivalent to

where $a, b, c$ and $d$ are the coefficients of the next-generation matrix (3.1). Next, note that (3.6) implies

Combining (5.21) with (5.20), we see that

Equation (5.20) says $(a-1)(d-1) > bc > 0$ and, therefore, $a-1$ and $d-1$ must have the same sign. Also, in view of (5.11), the condition $\mathscr{R}_0 < 1$ implies that $\frac{a + d} 2 < 1$, whereby $(a-1) + (d-1) < 0$ and it must be that $a-1 < 0$ or $d-1 < 0$. Having the same sign,

Finally, the first condition in (5.19) is $(a-1)\Upsilon_f + (d-1)\Upsilon_m < 0$, which follows immediately from (5.21). Thus, we conclude that the Disease Free Equilibrium of the model (2.3) is LAS if ${\mathscr{R}}_{0} < 1$ and unstable if ${\mathscr{R}}_0 > 1.$

We shall compute the non-trivial equilibria of our model assuming that there are no disease-induced deaths. An upper index $^*$ next to any state variable is used to denote an equilibrium value of that variable. It follows that $N_k(t)$ approaches the limiting population size $N_k^*$ given in (5.2). The theory of asymptotically autonomous systems, as discussed by Castillo-Chavez and Thieme ^[18], and Markus ^[19], indicates that if $N_k$ approaches a constant limit, then the system is asymptotically equivalent to one where $N_k$ is replaced by this limit.

To compute the nontrivial boundary equilibria ("only-one-sex-infected"), i.e., those where $I_f^* = 0$ and $I_m^* > 0$, or $I_m^* = 0$ and $I_f^* > 0$, we shall assume, without loss of generality, that $I_f^* = 0$ and $I_m^* > 0$. From the second equation in (2.3), we see that this condition is equivalent to $\beta_{mf} = 0$ and $\beta_{mm} > 0$.

Proposition 3.3. The "only-males-infected" equilibrium $E^*_m = \left(N^*_f, 0, 0, S^*_m, I^*_m, R^*_m\right)$, given by

exists in $\mathscr{D}$ and differs from the DFE if, and only if, ${\mathscr{R}}_{0, m} > 1$. It is locally asymptotically stable if, and only if, ${\mathscr{R}}_{0, m} > 1$ and ${\mathscr{R}}_{0, f} < 1$.

Proof. First, note that

From $\dot R_m = 0$ in (2.3), we obtain

Then, from $\dot I_m = 0$ in (2.3), we obtain

whereby

The last relation shows explicitly what fraction $\varepsilon$, $0\le\varepsilon < \frac{\gamma_m+\mu_m}{\gamma_m+\mu_m+\theta_m} < 1$, of the total male population is infected.

Note that when immunity is not permanent (i.e., $\theta_m > 0$), it is not possible for the whole population to be infected because $\varepsilon$ is bounded away from 1 no matter how large $\mathscr{R}_{0, m}$ might be.

Next, by (5.22)–(5.24), we have

Finally, by the third equation in (2.3) together with (5.24), we obtain

For the "Only-Females-Infected" Endemic Equilibrium, we can proceed analogously assuming $I_m^* = 0$ and $I_f^* > 0$ (equivalent to $\beta_{fm} = 0$ and $\beta_{ff} > 0$).

Proposition 3.4. The "only-females-infected" equilibrium $E^*_f = \left(S^*_f, I^*_f, R^*_f, N^*_m, 0, 0\right)$, given by

exists in $\mathscr{D}$ and differs from the DFE if, and only if, ${\mathscr{R}}_{0, f} > 1$. It is locally asymptotically stable if, and only if, ${\mathscr{R}}_{0, f} > 1$ and ${\mathscr{R}}_{0, m} < 1$.

Proof. The stability of $E_f^*$ can be studied analogously to that $E_m^*$, based on linearization of the model about $E^*_f$ and checking that all the eigenvalues of its Jacobian have negative real part if, and only if, ${\mathscr{R}}_{0, f} > 1$ and ${\mathscr{R}}_{0, m} < 1$.

Now, let us assume $I_f^*\not = 0$ and $I_m^*\not = 0$, and determine an algebraic system to find all endemic equilibria (i.e., in the interior of $\mathscr{D}$). To do this, we use the fact that the population of each sex $i = f, m$ is asymptotically constant of size $N^*_i$ and thus $S^*_i = N^*_i-(I^*_i+R^*_i).$

From $\dot R_f = 0$ in (2.3) we get

Moreover, from $\dot I_f^* = 0$ in (2.3) and (5.27) we have

where,

We next proceed similarly for male compartments. We have $\dot R_m = 0$ in (2.3), whereby

Moreover, from $\dot I_m^* = 0$ in (2.3) and (5.30) we have

where,

Thus, the infected classes at steady state must satisfy the following coupled two-by-two quadratic system:

We can find explicitly endemic equilibria when intrasex transmission occurs in only one direction (i.e., $\beta_{mf} = 0$ or $\beta_{fm} = 0$, but not both).

Proposition 3.5. Let $\mathscr{R}_{0, f} > 1$ and $\mathscr{R}_{0, m} > 1$ —whereby a fortiori $\mathscr{R}_0 > 1$ by Proposition 3.1. If $\beta_{fm} = 0$ (i.e., females cannot infect males), then the model admits a unique endemic equilibrium $E^* = (S^*_f, I^*_f, R^*_f, S^*_m, I^*_m, R^*_m)$, given by

Proof. From (5.28) and (5.29), we see that $b_2 = 0\iff\beta_{fm} = 0\iff b_4 = 0$, and then it must be either $I^*_f = 0$ or $I^*_f = \frac{b_1}{b_3}$, giving the size of the infected females compartment as that in the disease-free equilibrium (if $\mathscr{R}_{0, f} < 1$) or in the "only-females-infected" equilibrium (3.8) (if $\mathscr{R}_{0, f} > 1$). Then, from (5.31) we obtain the nontrivial size of the infected males compartment in this case as the only positive root of the quadratic equation $-a_3I^2+\left(a_1-{a_4I_f}\right)I+{a_2I_f} = 0$, namely $I_m^* = \frac{a_1}{a_3}$ when $I_f^* = 0$ (which is the case of "only-males-infected", $\mathscr{R}_{0, f} < 1$ and $\mathscr{R}_{0, m} > 1$) or, when $I^*_f = \frac{b_1}{b_3}$,

which is positive if $\mathscr{R}_{0, m} > 1$, by (5.32).

Proposition 3.6. Let $\mathscr{R}_{0, f} > 1$ and $\mathscr{R}_{0, m} > 1$ —whereby a fortiori $\mathscr{R}_0 > 1$ by Proposition 3.1. If $\beta_{mf} = 0$ (i.e., males cannot infect females), then the model admits a unique endemic equilibrium $E^* = (S^*_f, I^*_f, R^*_f, S^*_m, I^*_m, R^*_m)$, given by

Proof. From (5.31) and (5.32), we see that $a_2 = 0\iff\beta_{mf} = 0\iff a_4 = 0$, and then it must be either $I^*_m = 0$ or $I^*_m = \frac{a_1}{a_3}$, giving the size of the infected males compartment in this case as that in the disease-free equilibrium (if $\mathscr{R}_{0, m} < 1$) or in the "only-males-infected" equilibrium (3.7) (if $\mathscr{R}_{0, m} > 1$). Then, from (5.28) we obtain the nontrivial size of the infected females compartment as the only positive root of the quadratic $-b_3I^2+\left(b_1-{b_4I_m}\right)I+{b_2I_m} = 0$, namely $I_f^* = \frac{b_1}{b_3}$ when $I_m^* = 0$ (which is the case of "only-females-infected", $\mathscr{R}_{0, f} > 1$ and $\mathscr{R}_{0, m} < 1$) or, when $I^*_m = \frac{a_1}{a_3}$,

which is positive if $\mathscr{R}_{0, f} > 1$, by (5.29).

For the general case when all transmission pathways are available, we obtain from the first equation in (5.33)

Thus, for $\mathscr{R}_{0, m} > 1$, we have $I_f > 0$ if

Proposition 3.7. Endemic equilibria for model (2.3) exist in the following cases:

1) The boundary equilibrium $E^*_{m} = (S_m^*, I_m^*, R_m^*, S_f^*, 0, 0)$ exists and is locally asymptotically stable if $\mathscr{R}_{0, m} > 1$ and $\mathscr{R}_{0, f} < 1$.

2) The boundary equilibrium $E^*_{f} = (S_m^*, 0, 0, S_f^*, I_f^*, R_f^*)$ exists and is locally asymptotically stable if $\mathscr{R}_{0, m} < 1$ and $\mathscr{R}_{0, f} > 1$.

3) At least one positive equilibrium $E^* = (S_m^*, I_m^*, R_m^*, S_f^*, I_f^*, R_f^*)$ exists if at least one of $\sigma_0, \sigma_1$, or $\sigma_2$ given by (5.42) is bigger than 1, and at least one of the positive roots $I^*_m$ of (5.43) satisfies the bounds given in (5.37).

Proof. Cases 1 and 2 are Propositions 3.5 and 3.6. As for case 3, substituting (5.36) into the second equation of (5.33), we have

or, equivalently,

where

Here,

Therefore, in order to have a positive equilibrium, it is enough to find $I_m\in \mathscr{I}$, such that

According to Descartes' Rule of Signs (see, for example, ^[8]), the maximum number of positive real roots a polynomial may have is given by the number of sign changes in consecutive pairs of coefficients of the polynomial. Since $I^*_m$ is the root of a cubic, we may have up to three positive real roots. Because (5.41) implies $\alpha_3 < 0$, the only way for this to occur is that $\alpha_0 > 0$, $\alpha_1 < 0$, and $\alpha_2 > 0$. Moreover, the only way for no positive real roots to exist is that $\alpha_0, \alpha_1, \alpha_2 < 0$, that is, $\sigma_0, \sigma_1, \sigma_2 < 1$. Thus, by examining all possible sign combinations, we get the results shown in Table 3.