A shiny app for modeling the lifetime in primary breast cancer patients through phase-type distributions

1.   Introduction

The World Health Organization defines sexually transmitted diseases (STDs) as various diseases that are transmitted through sexual contact, similar sexual behaviors and indirect contact. The common reasons of STDs are bacteria, yeast and viruses^[1]. STDs such as Trichomoniasis, Gonorrhea, Syphilis, Genital Warts and Herpes have become a serious public health problem. Recently, mathematical models of epidemic or population dynamics have been widely used ^{[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]}.

The mathematical models in the early researches of STDs generally assume that both the males and females are evenly mixed, that is, the contacts of all individuals are equal. These models ignore the social and contact structures of the real population. For instance, the number of sexual partners in different individuals may vary^[21]. One method described in ^[22] that includes contact heterogeneity is the core group model. Although this model divides the population into two categories with a large number of sexual contacts and less sexual contact, it is still considered that the individuals are well mixed. Thus this model is not suitable for the spread of disease among the general public and is more suitable for sex workers.

People gradually have realized the importance of heterogeneous social networks in recent years. The spread of sexually transmitted diseases occurs in social networks based on real human contact. The so-called network contains many nodes representing different individuals in the real system and the edges of the connected nodes representing relationships between individuals. Nodes are often separated into two categories by sexual contacts, and only nodes of the opposite type can be connected. The contacts between the opposite sex are represented by a bipartite network^[23].

Most of the previous models on complex networks assume that disease transmission is a Poisson process, and every individual randomly selects an individual from the population. This assumption implies that the duration of the partnership is very short. The focus of many researches gradually begins to understand the role of some individuals those with many connections in two ways. One is assuming a short-lived partnership (the time of disease transmission is much longer than the duration of the partnership)^[24], another is assuming that the network is static (the time of disease transmission is much shorter than the duration of the partnership)^{[25,26,27,28]}. The edge-based compartmental model (EBCM) has the potential to unify these two approaches recently, and allows partnership durations to last from zero to infinity^[29,30,31].

The above models assumed that the initial infection ratio was infinitesimally small. The inapplicability of this assumption was $R_{0} < 1$ or the initial infection rate was not negligible^[25]. Miller^[32] extended the edge-based compartmental model to arbitrary initial conditions and gave a detailed explanation of the part of the initial proportion infected could not be ignored. This helps to resolve an obvious paradox in early work, that is, if there are too many people initially infected, the number of susceptible people may increase. This also helps to explain a significant small deviation observed between the simulation and theory in the previous paper^[25]. This modification makes sense for us to consider vaccination or previous infections. Yan et al.^[33] considered the spread of STDs SIR sexually transmitted diseases on bipartite networks representing heterosexual individuals.

Motivated by ^[32,33], some sexually transmitted diseases have the latent period but the individuals are less infectivity during the latent period, so we introduce latent compartment in our model and assume that transmission rate of the latent period is zero in this paper. We first assume that the proportion of initial infections is infinitesimally small. Based on this, the qualitative and stability of our model is further considered when the proportion of initial infections is arbitrarily large.

This paper consists of 6 parts. We derive the edge-based SEIR model for sexually transmitted diseases in section 2. We computer the reproduction number $R_{0}$ and the final epidemic size of the disease and analyze the local dynamics of the model in consideration of the infinitesimal initial infection rates in section 3. In section 4, we further analyze the dynamical behavior of our model considering a large number of infected individuals at the initial moment. In section 5, we perform some simulations with different initial values on different networks and some sensitivity analysis. The final section of the paper gives some concluding remarks.

2.   Model derivation

In the section, the network during the epidemic is assumed to be fixed and be of configuration type^[34], and disease deaths are ignored. In the network, $(P_M(k))$ and $(P_F(k))$ represent the distribution of male and female individuals, respectively. Following^[29,33], We distribute stubs to every men and women according to $(P_M(k))$ and $(P_F(k))$ at random. Further we pick out two stubs attached to individuals of different genders and connect them. We keep repeating the process until no new edge appear. Multiple edges, degree correlation, self loops, and clustering are negligible on the network constructed by this method^[33,35].

2.1.   Variables and parameters

In Table 1, we list some variables and parameters. $U_{M} (U_{F})$ is a male (female) individual being tested which is randomly selected at the initial moment. The proportion of individuals in a certain state in the population is equal to the probability that $U_{M} (U_{F})$ in this state. We modify $U_{M} (U_{F})$ so that he does not be transmitted to any of its partners when infected. More discussion about the tested male $U_{M}$ and female $U_{F}$ is in ^[29,36]. $S_M, E_M, I_M, R_M$ and $S_F, E_F, I_F, R_F$ are the proportion of the susceptible, exposed, infected, recovered individuals in male and female individuals, respectively. They also are the probabilities of the tested male $U_{M}$ and the tested female $U_{F}$ are susceptible, exposed, infected or recovered, respectively. Other variables and parameters can be found in Table 1.

2.2.   Equation derivation

In this section, we assume that individuals in exposed compartment are not infectivity, then deduce a system portraying the spread of sexually transmitted diseases including above variables and parameters. The fraction of the male individuals who are susceptible at time t $S_M(t)$ is derived at first. $S_M(k, 0)$ is the probability which the tested male $U_M$ is susceptible and has degree k at t = 0, and $\theta_{M}(t)$ is the probability that he has not been infected by his partner for a period of time. So we have

So, the proportion of the male individuals who are susceptible at time t is

We can get $S_M(t)$ and $S_F(t)$ by the equations of $\theta_M$ and $\theta_F$.

We can obtain $\dot{I}_M = \upsilon_M(1-S_M-I_M-R_M)-\gamma_MI_M$ by combining equations $\dot{I}_M = \upsilon_ME_M-\gamma_MI_M$ and $S_M+E_M+I_M+R_M = 1$. We also know that $R_M$ satisfy $\dot{R}_M = \gamma_MI_M$ and $S_M(t) = \Psi_M(\theta_M(t))$, so we can completely define $S_M, E_M, I_M$ and $R_M$ assuming $\theta_M(t)$ and initial conditions for $R_M$ and $I_M$ are known. Similarly, we can completely define $S_F, E_F, I_F$ and $R_F$ assuming $\theta_F(t)$ and initial conditions for $R_F$ and $I_F$ are known.

If we get the equations of $\theta_M$ and $\theta_F$, we can close the system. We already know $\theta_M(t)$ is the probability that a tested male individual has not been infected by his randomly chosen partner yet. These partners are made up of susceptible, exposed, infected and recovered females, so we have $\theta_M = \phi_{S_{F}M}+\phi_{E_{F}M}+\phi_{I_{F}M}+\phi_{R_{F}M}.$ Because $\phi_{I_{F}M}$ is the probability of a randomly chosen partner of $U_M$ who does not transmit the disease to $U_M$ before is infectious at time t, we have

Similarly, we have

Now we need to get the equations of $\phi_{I_{F}M}$ and $\phi_{I_{M}F}$. Because we have the equation $\theta_M = \phi_{S_{F}M}+\phi_{E_{F}M}+\phi_{I_{F}M}+\phi_{R_{F}M}.$ We can find the $\phi_{S_{F}M}$ class, noticing the probability that $U_M$ has a female partner who is susceptible at the initial moment is $\phi_{S_{F}M}(0)$ and the probability of the susceptible female has degree k is $\frac{kP_{F}(k)S_F(k, 0)}{\sum_{j}jP_{F}(j)S_F(j, 0)}$. So her probability of being a susceptible individual after time t is $\sum_k\frac{kP_{F}(k)S_F(k, 0)\theta_F^{k-1}}{\sum_{j}jP_{F}(j)S_F(j, 0)}$. Thus we have

Next we start to calculate the $\phi_{R_{F}M}$ class, From (Figure 1) we notice that only one edge enters the $\phi_{R_{F}M}$ class. We have

Integrating Eq (2.1) and Eq (2.2), we have

Integrating Eq (2.3) from 0 to t yields

So we have

In addition, in the Figure 1 we also notice that there are two edges leaving the class $\phi_{I_{F}M}$ at rates $\beta_{FM}$ and $\gamma_{F}$, respectively. And an edge enters the class $\phi_{I_{F}M}$. We have

At last, the following equation can be obtained,

In summary, our model can be derived to be

We can further simplify the model by substituting Eq (2.4) into Eq (2.5), we have

Similarly, we have

Now we can rewrite the model (2.6) as

3.   Analysis of the model

3.1.   Disease free equilibrium and the basic reproduction number

Considering the following equations of model (2.6) :

It is easy to know that $(\phi_{E_{F}M}, \phi_{E_{M}F}, \phi_{I_{F}M}, \phi_{I_{M}F}, \theta_M, \theta_F) = (0, 0, 0, 0, 1, 1)$ is the disease free equilibrium of system (3.1). We calculate the basic reproduction number $R_0$ by applying the method of the second generation matrix in ^[37]. In our model (3.1), the classes $\phi_{E_{F}M}$ and $\phi_{E_{M}F}$ act as "exposed" types and the classes $\phi_{I_{F}M}$ and $\phi_{I_{M}F}$ act as "infected" types. Variables $\theta_{M}$ and $\theta_{F}$ act as "susceptible" types since they can enter the $\phi_{E_{F}M}$ and $\phi_{E_{M}F}$ classes when the disease breaks out. So we only need to linearize these equations about $\phi_{E_{F}M}$, $\phi_{E_{M}F}$, $\phi_{I_{F}M}$ and $\phi_{I_{M}F}$ in (3.1) at the disease free equilibrium $(\phi_{E_{F}M} = \phi_{E_{M}F} = \phi_{I_{F}M} = \phi_{I_{M}F} = 0)$, we have

Applying the method of the second generation matrix, we obtain

where

and

Thus,

Hence,

where $\rho$ represents the spectral radius and the $R_{0}$ is the basic reproduction number. The biological interpretation of $\frac{\Psi''_{M}(1)}{\Psi'_{M}(1)}$ comes from observing the situation of a random individual's partner $V_M$ in the early stages of the epidemic. If $V_M$ is infected by that randomly infected female, then $\frac{\Psi''_{M}(1)}{\Psi'_{M}(1)}$ is the expectant number of other partners $V_M$ has (his excess degree). Then $\frac{\beta_{MF}}{{\beta_{MF}}+\gamma_M}$ is the possibility that an infected male individual transmit the disease to his partner. So we can get $\frac{\phi_{S_{M}F}(0)\beta_{MF}\Psi''_{M}(1)}{(\beta_{MF}+\gamma_M)\Psi'_{M}(1)}$ is the number of individuals who may be infected by $V_M$. We have a similar result for $V_F$ infected with male individual who are randomly infected, that is, $\frac{\phi_{S_{F}M}(0)\beta_{FM}\Psi''_{F}(1)}{(\beta_{FM}+\gamma_F)\Psi'_{F}(1)}$ is the number of individuals who may be infected by $V_F$. So $R_0$ is the geometric mean of the number of individuals infected with $V_M$ and the number of individuals infected with $V_F$, which is consistent with the result we calculated.

We usually calculate basic reproduction number of the disease at infinitesimal initial values $(\phi_{S_{F}M}(0) = \phi_{S_{M}F}(0) = 1$, $\phi_{R_{F}M}(0) = \phi_{R_{M}F}(0) = 0)$, then $R_{0}$ in Eq (3.2) becomes

Yan et al.^[33] derived the basic reproduction number of a new edge-based SIR model of sexually transmitted diseases on bipartite networks. Comparing to that of ^[33], we know that the $R_0$ of our model is the same. We also note that $R_0$ is symmetric in the parameters describing of male and female properties being consistent with ^[38].

3.2.   The final epidemic size

The final size relation has been done for various models with small or large initial conditions ^[38,36]. Motivated by ^[38,36], the final epidemic size of our model be derived in what follows.

We set $\frac{\mbox{d}}{\mbox{d}t}\theta_{M} = \frac{\mbox{d}}{\mbox{d}t}\theta_{F} = \frac{\mbox{d}}{\mbox{d}t}\phi_{I_{F}M} = \frac{\mbox{d}}{\mbox{d}t}\phi_{I_{M}F} = \frac{\mbox{d}}{\mbox{d}t}I_{F} = \frac{\mbox{d}}{\mbox{d}t}I_{M} = \frac{\mbox{d}}{\mbox{d}t}E_{F} = \frac{\mbox{d}}{\mbox{d}t}E_{M} = 0$, so $\phi_{E_{F}M}(\infty) = \phi_{E_{M}F}(\infty) = E_F(\infty) = E_M(\infty) = 0$ and $\phi_{I_{F}M}(\infty) = \phi_{I_{M}F}(\infty) = I_F(\infty) = I_M(\infty) = 0$. From (2.7), we have

and

Since $\frac{\Psi'_{F}(\theta_{F})}{\Psi'_{F}(1)} = \sum_k\frac{kP_{F}(k)S_F(k, 0)\theta_F^{k-1}}{\sum_{j}jP_{F}(j)S_F(j, 0)}$ and $\frac{\Psi'_{M}(\theta_{M})}{\Psi'_{M}(1)} = \sum_k\frac{kP_{M}(k)S_M(k, 0)\theta_M^{k-1}}{\sum_{j}jP_{M}(j)S_M(j, 0)}$, we have

and

Then

and

Since we have $S_M(\infty) = \Psi_M(\theta_M(\infty))$ and $S_F(\infty) = \Psi_F(\theta_F(\infty))$, we can get the final epidemic size with arbitrary initial conditions are

and

We calculate the final epidemic size of the disease at infinitesimal initial values, that is, $\phi_{R_{F}M}(0) = \phi_{R_{M}F}(0) = 0, \phi_{S_{F}M}(0) = \phi_{S_{M}F}(0) = 1, R_M(0) = R_F(0) = 0$. Then,

and

Further we get the final epidemic size at infinitesimal initial values

and

From Eqs (3.2), (3.3) and (3.4), we know that the basic reproduction number and the final size of an epidemic are not related to the length of latent period in the SEIR model without infectivity during the latent period.

3.3.   Stability of disease free equilibrium with infinitesimal initial conditions

In the section, we investigate the disease equilibrium of our model at the infinitesimal initial conditions, i.e. $\phi_{E_{F}M}(0) = \phi_{I_{F}M}(0) = 0$ and $\phi_{E_{M}F}(0) = \phi_{I_{M}F}(0) = 0$ or $\phi_{S_{F}M}(0)+\phi_{R_{F}M}(0) = 1$ and $\phi_{R_{M}F}(0)+\phi_{S_{M}F}(0) = 1$. We only need to study the equation group consisting of the equations of $\phi_{I_{F}M}, \phi_{I_{M}F}, \theta_M$ and $\theta_F$ in the model (2.7), i.e.

Noting that $\phi_{I_{F}M}, \phi_{I_{M}F}, \theta_M$ and $\theta_F$ are all probabilities, we only need to consider this system in $\Omega = \{(\phi_{I_{F}M}, \phi_{I_{M}F}, \theta_M, \theta_F)| 0\leq\phi_{I_{F}M}, \phi_{I_{M}F}, \theta_M, \theta_F\leq1\}$. It is easy to verify that $\Omega$ is a positive invariant set of system (3.5). We have these results in what follows.

Theorem 1. There is a disease free equilibrium $E^0(0, 0, 1, 1)$ in the system (3.5) with infinitesimal infected initial values. Moreover,

(Ⅰ) if $R_0 < 1$, the disease free equilibrium $E^0$ is locally asymptotically stable,

(Ⅱ) if $R_0 > 1$, there exists only one endemic equilibrium $E^* = (0, 0, \theta_M^{*}, \theta_F^{*})$ in which $0 < \theta_M^{*}, \theta_F^{*} < 1$, and it is locally asymptotically stable.

Proof. We know that the equilibria need to satisfy the following.

Now we study system (3.5) in $\Omega$. It is easy to know that $E^0 = (0, 0, 1, 1)$ is always a disease free equilibrium in the system (3.5). Below we use three steps to complete our proof.

(1) Noting that for system (3.5) whose initial infection is infinitesimal, the characteristic equation at the disease free equilibrium $E^0$ is

Therefore, Eq (3.6) can be written as

where

Applying Routh-Hurwitz criteria, if $R_0 < 1$, we have $\frac{\beta_{FM}\beta_{MF}\phi_{S_{M}F}(0)\phi_{S_{F}M}(0)\Psi''_{M}(1)\Psi''_{F}(1)}{(\beta_{FM}+\gamma_F)(\beta_{MF}+\gamma_M)\Psi'_{M}(1)\Psi'_{F}(1)} < 1,$ further we can get

So $a_{4} > 0$, and we have

For the sake of clarity, we put the formulae of $b_{i}$($i = 1, 2, 3, 4$) in the Appendix A. We can clearly know that the disease free equilibrium $E^0$ is locally asymptotically stable if $R_0 < 1$.

(2) We show that there is a unique endemic equilibrium $E^* = (0, 0, \theta_M^{*}, \theta_F^{*})$ in which $0 < \theta_M^{*}, \theta_F^{*} < 1$ when $R_0 > 1$. Inspired by ^[33], we know that $\phi_{I_{F}M}$ and $\phi_{I_{M}F}$ are always equal to 0 at the equilibrium and construct the auxiliary function from system (2.7):

where

We only need to prove that $f(\theta_M)$ has a unique solution between 0 and 1. It is easy for us to compute that $f(0) < 0$ and $f(1) = 0$, we have

and

According to $\frac{d^{2}f(\theta_M)}{\mbox{d}\theta_M^2}\leq0$, we find that $f(\theta_M)$ is concave. Therefore $\frac{\mbox{d}f(\theta_M)}{\mbox{d}\theta_M}$ is monotonically decreasing in the interval (0, 1). We have $\frac{\beta_{FM}\phi_{S_{F}M}(0)}{\beta_{FM}+\gamma_F}\frac{\Psi''_{F}(1)}{\Psi'_{F}(1)}\frac{\beta_{MF}\phi_{S_{M}F}(0)}{\beta_{MF}+\gamma_M}\frac{\Psi''_{M}(1)}{\Psi'_{M}(1)} > 1$ according $R_{0} > 1$, we also find $\frac{\mbox{d}f(0)}{\mbox{d}\theta_M} = 1 > 0$ and $\frac{\mbox{d}f(1)}{\mbox{d}\theta_M} < 0$ from Eq (3.7). Therefore, we can derive a $\theta_M^{**}$ in the interval (0, 1) when $R_{0} > 1$ so that $\frac{\mbox{d}f(\theta_M^{**})}{\mbox{d}\theta_M} = 0$ according to the intermediate value theorem. Since $f(\theta_M)$ is a concave function in the interval (0, 1), we can get $\theta_M^{**}$ to make $f(\theta_M)$ the largest in the interval, so $\theta_M^{**}\in(0, 1)$ and $f(\theta_M^{**}) > f(1) = 0$. Considering the monotonicity of $f(\theta_M)$, we use the intermediate value theorem to find a $\theta_M^{*}\in(0, \theta_M^{**})$ makes $f(\theta_M^{*}) = 0$ and $\theta_M^{*}$ is unique on the basis of the increment of $f(\theta_M)$ in $(0, \theta_M^{**})$ and the decrement of $f(\theta_M)$ in $(\theta_M^{**}, 1)$ (see Figure 2).

(3) We prove that the endemic equilibrium $E^*$ is locally asymptotically stable in $\Omega$. The characteristic equation at the endemic equilibrium $E^*$ is

Therefore, Eq (3.9) can be written as

where

Applying Routh-Hurwitz, we have

For the sake of clarity, we put formulae of $d_{i}$ ($i = 1, 2, 3, 4$) in the Appendix B. We can clearly know that the endemic equilibrium $E^*$ is locally asymptotically stable. The proof is completed.

4.   The model with large initial conditions

We analyze the local dynamics of system (2.7) with larger initial value of infection in this section. The null lines of $\theta_{M}$ and $\theta_{F}$ in the model (2.7) are in what follows:

Where $\phi_{R_{M}F}(0)$ and $\phi_{R_{F}M}(0)$ are fixed, as the values of $\phi_{I_{M}F}(0)$ and $\phi_{I_{F}M}(0)$ increase from 0, $L_{M}$ and $L_{F}$ move to the left and down, respectively. We have the following two claims:

(1) If $R_0 < 1$, system (2.7) for infinitesimal initial conditions have only one locally asymptotically stable disease free equilibrium $E^0(0, 0, 1, 1)$ in $\Omega$. As $\phi_{I_{M}F}(0)$ and $\phi_{I_{F}M}(0)$ gradually increase from 0, the disease free equilibrium $E^0$ moves from (0, 0, 1, 1) to the lower left to an internal point $E^0_{*} = (0, 0, \theta_{M*}^{(0)}, \theta_{F*}^{(0)})$ of $\Omega$ in which $0 < \theta_{M*}^{(0)}, \theta_{F*}^{(0)} < 1$, and $E^0_{*}$ and $E^0$ are consistent in local stability.

(2) If $R_0 > 1$, system (2.7) for infinitesimal initial conditions has one locally asymptotically stable endemic equilibrium $E^* = (0, 0, \theta_M^{*}, \theta_F^{*})$ and an unstable disease free equilibrium $E^0(0, 0, 1, 1)$ in $\Omega$. As $\phi_{I_{M}F}(0)$ and $\phi_{I_{F}M}(0)$ gradually increase from 0, disease free equilibrium $E^0$ moves from (0, 0, 1, 1) to the upper right to an external point $E^0_{*} = (0, 0, \theta_{M*}^{(0)}, \theta_{F*}^{(0)})$ with min$\{\theta_{M*}^{(0)}, \theta_{F*}^{(0)}\} > 1$ (we do not analyze its dynamics since this point is not in $\Omega$), and the endemic equilibrium $E^*$ moves from $(0, 0, \theta_M^{*}, \theta_F^{*})$ slightly to the lower left to another internal point $E^*_{*} = (0, 0, \theta_{M*}^{*}, \theta_{F*}^{*})$ of $\Omega$ in which $0 < \theta_{M*}^{(*)} < \theta_M^{*}$ and $0 < \theta_{F*}^{(*)} < \theta_F^{*}$. In addition. The stability of $E^0_{*}$ and $E^*_{*}$ is consistent with $E^0$ and $E^*$ respectively.

In summary, we give the following theorem for system (2.7) with max$\{\phi_{I_{M}F}(0), \phi_{I_{F}M}(0)\} > 0$.

Theorem 2. For our system (2.7) with large initial values of infection.

(Ⅰ) if $R_0 < 1$, system has only one disease free equilibrium $E^0_{*} = (0, 0, \theta_{M*}^{(0)}, \theta_{F*}^{(0)})$ of $\Omega$ with $0 < \theta_{M*}^{(0)}, \theta_{F*}^{(0)} < 1$, and the solution gradually approaches $(0, 0, \theta_{M*}^{(0)}, \theta_{F*}^{(0)})$ from (0, 0, 1, 1),

(Ⅱ) if $R_0 > 1$, system has only one endemic equilibrium $E^*_{*} = (0, 0, \theta_{M*}^{*}, \theta_{F*}^{*})$ in $\Omega$, and the solution gradually approaches $(0, 0, \theta_{M*}^{*}, \theta_{F*}^{*})$ from (0, 0, 1, 1).

Proof. Similar to the proof of Theorem 1, we show that there is a unique equilibrium in $\Omega$. Since $\phi_{I_{F}M}$ and $\phi_{I_{M}F}$ are always equal to zero at the equilibrium point, we construct the following auxiliary function:

where

We only need to prove that there is a unique solution for $g(\theta_M)$ between 0 and 1. Getting $g(0) < 0$ is easy for us, and we have

and

According to $\frac{d^{2}g(\theta_M)}{\mbox{d}\theta_M^2}\leq0$ for all $\theta_M\geq0$ in Eq (4.1), we find that $g(\theta_M)$ is concave. When the initial infection is arbitrary large, we have

$0 < \phi_{S_{F}M}(0)+\phi_{R_{F}M}(0) < 1, 0 < \phi_{S_{M}F}(0)+\phi_{R_{M}F}(0) < 1, \theta_M(0) = 1, \theta_F(0) = 1.$

According to $0 < \theta_F < 1$, we have

Considering $g(\theta_M)$ is concave when $\theta_M\geq0$, we use the intermediate value theorem to find that there is a unique $\tilde{\theta_M}$ in (0, 1) to make $g(\tilde{\theta_M}) = 0$. And we have

We conclude that our system (2.7) has only one equilibrium $\tilde{E} = (0, 0, \tilde{\theta_M}, \tilde{\theta_F})$ in $\Omega$. It is easy to conclude that $\tilde{E}$ is locally asymptotically stable which is similar to the proof process of Theorem 1. And the solutions of our system gradually approach $\tilde{E}$ from (0, 0, 1, 1). The proof is completed.

Example 1. We assume that both males and females in system (2.7) obey the Poisson distributions given by $P_M(k) = P_F(k) = \frac{\lambda^{k}e^{-\lambda}}{k!}$ with $\lambda = 5,$ $\phi_{R_{F}M}(0) = \phi_{R_{M}F}(0) = 0.$ We get $R_0 = 0.5556 < 1$ and a unique locally asymptotically stable disease free equilibrium $E^{0} = (0, 0, 1, 1)$ in $\Omega$ (see Figure 3(a)). Moreover, we also get $R_0 = 0.25 < 1$ and a unique locally asymptotically stable disease free equilibrium $E^0_{*} = (0, 0, 0.9070, 0.9071)$ in $\Omega$ (see Figure 4(a)). It is obvious that the solution gradually approaches $(0, 0, 0.9070, 0.9071)$ from (0, 0, 1, 1) as $\phi_{I_{F}M}(0)$ and $\phi_{I_{M}F}(0)$ increase from 0.

We set parameter values and initial values, getting $R_0 = 2 > 1$, a disease free equilibrium $E^0$ and the locally asymptotically stable endemic equilibrium $E^{*} = (0, 0, 0.6813, 0.6813)$ in $\Omega$ (see Figure 3(b). Moreover, we also get $R_0 = 1.8 > 1$ and a unique locally asymptotically stable endemic equilibrium $E^*_{*} = (0, 0, 0.6763, 0.6765)$ in $\Omega$ (see Figure 4(b)). It is obvious that the solution gradually approaches $(0, 0, 0.6763, 0.6765)$ from (0, 0, 1, 1) as $\phi_{I_{F}M}(0)$ and $\phi_{I_{M}F}(0)$ increase from 0.

5.   Stochastic simulation and sensitivity analysis

We implement the comparison between stochastic simulations and numerical predictions to our SEIR model on Poisson and scale-free networks in this section. The distributions of Poisson network and scale-free network are given by $P(k) = \frac{\lambda^{k}e^{-\lambda}}{k!}(1 \leq k \leq10)$ and $P(k) = (r-1)m^{(r-1)}k^{-r}(3 \leq k \leq10)$, respectively. The average degree of these two networks is 5. For network degree distributions of male and female individuals, we use the configuration model described in Section 2 to generate random contact.

Our model is defined in a random simulation as follows. The model has no infectivity during the latent period. Once a susceptible male is infected with a rate of $\beta_{FM}i_F$, where $i_F$ is the partner (female) node which he is exposed to, and he first passes through the latent period with an average length of $1/\upsilon_M$, then enters the infectious class. The same is true for female nodes. We define parameters are $\beta_{FM} = \beta_{MF} = 0.01, \gamma_F = \gamma_M = 0.04, \upsilon_F = \upsilon_M = 0.2$. The initial values are $\phi_{E_{F}M}(0) = \phi_{E_{M}F}(0) = 0.01, \phi_{I_{F}M}(0) = \phi_{I_{M}F}(0) = 0.04, \phi_{R_{F}M}(0) = \phi_{R_{M}F}(0) = 0.06, E_M(0) = E_F(0) = 500, I_M(0) = I_F(0) = 200, R_M(0) = R_F(0) = 300.$ If there are no exposed and infected individuals in the networks, the entire spread of epidemic will stop. As we can see from Figures 5 and 6, the prediction of model and average of stochastic simulation fit well on both types of networks with the same average degree. This shows that model we built can accurately simulate the spread of disease.

Considering that many sexually transmitted diseases have different proportions of exposed and infectious individuals among male and female populations. We set the initial values of the male population different from the female population on two networks in Figures 7 and 8. We observe that when the initial infection of women is lower than that of men, the peak of female infectious individuals will be higher than that of male. This result is consistent on both Poisson and scale-free networks. We can find that this is consistent with the real situation. When the number of female individuals infected is small initially, the number of male individuals infected by female will be smaller, and finally the peak of male individuals infected is lower than that of female.

In Figures 9-11, we change the initial conditions on the Poisson and the scale-free networks, and observe the dynamics of the model. In Figure 9(a), we change the value of the initial infectious individuals $I_M(0)( = I_F(0))$ with $R_M(0)( = R_F(0))$ is fixed, observing that the initial values change does not affect the peak arrival time of disease, but it affects the peak size of disease on different networks (i.e. the larger $I_M(0)( = I_F(0))$, the greater the peak of the disease). In Figure 9(b), we change the value of the initial recoverers $R_M(0)( = R_F(0))$ with $I_M(0)( = I_F(0))$ is fixed, observing that the change in the initial values does not affect the disease peak arrival time, but it affects the peak size on different networks (i.e. the larger $R_M(0)( = R_F(0))$, the greater the peak of the disease).

Figure 10(a) and (b) are contour maps of the initial infectious individuals on Poisson and scale-free networks respectively, then we observe that the change in $I_M(0)( = I_F(0))$ does not affect the peak arrival time of the infection but affects the peak value with that peak size is proportional to $I_M(0)$ and $I_F(0)$ on two types of networks. In addition, we also find that the scale-free network has a larger disease peak than the Poisson network in the same initial infections.

Figure 11(a) and (b) are contour maps of the initial recoverers on Poisson and scale-free networks respectively, then we find that the initial recoverers $R_M(0)( = R_F(0))$ affect the peak arrival time on two types of networks. It is obvious that the larger the $R_M(0)$ and $R_F(0)$ are, the earlier the peak arrives. We also find that the scale-free network has a larger disease peak than the Poisson network in the same initial recoverers.

From Figure 12, we know that changing the length of the latent period $1/\upsilon_M( = 1/\upsilon_F)$ can affects both the peak size of the infection and its arrival time. The shorter the latent period is, the larger the peak value an the earlier the arrival time. So we can take some measures to regulate the length of the latent period to interfere with the spread of the disease.

Figure 13(a) shows the final epidemic size $R_{M}(\infty)$. Parameters $\rho$ and $\langle k \rangle$ represent the ratio of the initial infected and exposed individuals and the average degree of the males and females respectively. We assume that the degree distribution in system (2.7) is subject to the Poisson network ($P(k) = \frac{\lambda^{k}e^{-\lambda}}{k!}(1 \leq k \leq20)$). We have $\Psi_M(x) = \Psi_F(x) = e^{\lambda(x-1)}$ and $\Psi'_{M}(1) = \Psi'_{F}(1) = \lambda = \langle k \rangle$. We set initial values $\phi_{R_{F}M}(0) = \phi_{R_{M}F}(0) = 0$ and parameters $\beta_{FM} = \beta_{MF} = 0.01, \gamma_F = \gamma_M = 0.04, \upsilon_M = \upsilon_F = 0.2$. We have $\phi_{S_{F}M}(0) = \phi_{S_{M}F}(0) = 1-\rho$. We can see from the figure that the final epidemic size increases when the average degree $\langle k \rangle$ increases, where $\rho$ remains unchanged. And when $\langle k \rangle$ remains unchanged, the final epidemic size increases with the increase of $\rho$. That is, the final epidemic size is proportional to $\rho$ and $\langle k \rangle$.

In Figure 13(b), we study the effect of each parameter on $R_0$. The PRCCs are calculated with respect to $\beta_{FM}, \beta_{MF}, \gamma_F, \gamma_M, \phi_{S_{F}M}(0), \phi_{S_{M}F}(0)$ and $\langle k \rangle$ with 2000 simulations. The input variables are subject to uniform distribution, and the positive and negative signs indicate that the effect is positive or negative, respectively. Sensitivity between 0 and 0.2 indicates that the parameter is weakly correlated, 0.2 to 0.4 is moderately correlated, and above 0.4 is highly correlated. As shown in Figure 13(b), we can see that $\beta_{FM}, \beta_{MF}, \phi_{S_{F}M}(0)$ and $\phi_{S_{M}F}(0)$ have a positive influence and are highly correlated on $R_0$, $\gamma_F$ and $\gamma_M$ have a negative influence and are highly correlated on $R_0$, and $\langle k \rangle$ is an insensitive parameter.

6.   Discussion and conclusion

In this paper, we extend an edge-based sexually transmitted SEIR model with no infectivity during the latent period, which the relationship between individuals is described by a bipartite network. We assume that the contact network is static and ignore the birth, death and migration of the population. We derive the basic reproduction number and the implicit formulas of the final epidemic size. We further analyze the dynamics of our model on different initial conditions, and observe the effects of different initial values on disease epidemics numerically. The basic reproduction number is consistent with that of ^[33] if we assume that transmission rate of the latent period is zero. Furthermore we also find that the length of the latent period is very much related to the arrival time and size of disease peak. How to construct and study the edge-based sexually transmitted models when the contact network is dynamic and the birth, death and migration of the population are considered are interesting. We leave this work in future.

Acknowledgments

We are grateful to the anonymous referees and the editors for their valuable comments and suggestions which improved the quality of the paper. This work is supported by the National Natural Science Foundation of China (11861044 and 11661050), and the HongLiu first-class disciplines Development Program of Lanzhou University of Technology.

Conflict of interest

The authors declare there is no conflict of interest.