Processing math: 59%
Research article Special Issues

Estimation of left ventricular parameters based on deep learning method


  • Estimating material properties of personalized human left ventricular (LV) modelling is a central problem in biomechanical studies. In this work we use deep learning (DL) method to evaluating the passive myocardial mechanical properties inversely. In the first part of the paper, we establish a standardized geometric model of the LV. The geometric model parameters are optimized based on 27 different healthy volunteers. In the second part, we use statistical methods and Latin hypercube sampling (LHS) to obtain the geometric parameters data. The LV myocardium is described using a structure-based orthotropic Holzapfel-Ogden constitutive law. The LV diastolic pressure-volume (PV) curves are calculated by numerical simulation. Tn the third part, we establish the multiple neural networks to pblackict PV curve parameters. Then, instead of using constrained optimization problems to solve constitutive parameters, DL was used to establish the nonlinear mapping relationship of geometric parameters, PV curve parameters and constitutive parameters. The results show that the deep learning method can greatly improve the computational efficiency of numerical simulation and increase the possibility of its application in rapid feedback of clinical data.

    Citation: Li Cai, Jie Jiao, Pengfei Ma, Wenxian Xie, Yongheng Wang. Estimation of left ventricular parameters based on deep learning method[J]. Mathematical Biosciences and Engineering, 2022, 19(7): 6638-6658. doi: 10.3934/mbe.2022312

    Related Papers:

    [1] Shubham Chaudhry, Gauri Agrawal, Maia Martcheva, A. K. Misra . Modeling the impact of temperature on the dynamics of carrier-dependent infectious diseases with control strategies. Mathematical Biosciences and Engineering, 2025, 22(7): 1722-1750. doi: 10.3934/mbe.2025063
    [2] Abdallah Alsammani, Calistus N. Ngonghala, Maia Martcheva . Impact of vaccination behavior on COVID-19 dynamics and economic outcomes. Mathematical Biosciences and Engineering, 2025, 22(9): 2300-2338. doi: 10.3934/mbe.2025084
    [3] Manoj Kumar, Syed Abbas, Abdessamad Tridane . Optimal control and stability analysis of an age-structured SEIRV model with imperfect vaccination. Mathematical Biosciences and Engineering, 2023, 20(8): 14438-14463. doi: 10.3934/mbe.2023646
    [4] Pride Duve, Samuel Charles, Justin Munyakazi, Renke Lühken, Peter Witbooi . A mathematical model for malaria disease dynamics with vaccination and infected immigrants. Mathematical Biosciences and Engineering, 2024, 21(1): 1082-1109. doi: 10.3934/mbe.2024045
    [5] Khadiza Akter Eme, Md Kamrujjaman, Muntasir Alam, Md Afsar Ali . Vaccination and combined optimal control measures for malaria prevention and spread mitigation. Mathematical Biosciences and Engineering, 2025, 22(8): 2039-2071. doi: 10.3934/mbe.2025075
    [6] Martin Luther Mann Manyombe, Joseph Mbang, Jean Lubuma, Berge Tsanou . Global dynamics of a vaccination model for infectious diseases with asymptomatic carriers. Mathematical Biosciences and Engineering, 2016, 13(4): 813-840. doi: 10.3934/mbe.2016019
    [7] Lan Zou, Jing Chen, Shigui Ruan . Modeling and analyzing the transmission dynamics of visceral leishmaniasis. Mathematical Biosciences and Engineering, 2017, 14(5&6): 1585-1604. doi: 10.3934/mbe.2017082
    [8] Jiazhe Lin, Rui Xu, Xiaohong Tian . Transmission dynamics of cholera with hyperinfectious and hypoinfectious vibrios: mathematical modelling and control strategies. Mathematical Biosciences and Engineering, 2019, 16(5): 4339-4358. doi: 10.3934/mbe.2019216
    [9] Shanjing Ren, Lingling Li . Global stability mathematical analysis for virus transmission model with latent age structure. Mathematical Biosciences and Engineering, 2022, 19(4): 3337-3349. doi: 10.3934/mbe.2022154
    [10] Lusha Shi, Jianghong Hu, Zhen Jin . Dynamics analysis of strangles with asymptomatic infected horses and long-term subclinical carriers. Mathematical Biosciences and Engineering, 2023, 20(10): 18386-18412. doi: 10.3934/mbe.2023817
  • Estimating material properties of personalized human left ventricular (LV) modelling is a central problem in biomechanical studies. In this work we use deep learning (DL) method to evaluating the passive myocardial mechanical properties inversely. In the first part of the paper, we establish a standardized geometric model of the LV. The geometric model parameters are optimized based on 27 different healthy volunteers. In the second part, we use statistical methods and Latin hypercube sampling (LHS) to obtain the geometric parameters data. The LV myocardium is described using a structure-based orthotropic Holzapfel-Ogden constitutive law. The LV diastolic pressure-volume (PV) curves are calculated by numerical simulation. Tn the third part, we establish the multiple neural networks to pblackict PV curve parameters. Then, instead of using constrained optimization problems to solve constitutive parameters, DL was used to establish the nonlinear mapping relationship of geometric parameters, PV curve parameters and constitutive parameters. The results show that the deep learning method can greatly improve the computational efficiency of numerical simulation and increase the possibility of its application in rapid feedback of clinical data.



    Uterine cervical cancer is a worldwide health problem but it is especially concerning in developing countries. It is the first or second most common cancer in women [1]. It is estimated that the probability of a person being infected with human papillomavirus (HPV) in their lifetime reaches 70 to 80% [2], and the total infection rate in the global population is as high as 11.7% [3]. An estimated 233,000 deaths were attributed to HPV infection in the year 2000 [4]. There were approximately 500,000 cases and 275,000 deaths due to cervical cancer worldwide in 2002, equivalent to about a tenth of all deaths in women due to cancer [5]. The burden of cervical cancer is disproportionately high (>80%) in the developing world [6].

    HPV was discovered to be the causative agent of cervical cancer in the 1970s by the Zur Hausen group [7]. Usually, the infecting papillomavirus is eliminated from individuals; however, some individuals retain the virus. Persistent infection with oncogenic HPV is recognized as the major cause of uterine cervical cancer [8]. Cervical carcinogenesis is a complex stepwise process over a continuum of increasingly severe precancerous changes known collectively as cervical intraepithelial neoplasia (CIN) [9]. The spectrum of CIN is traditionally divided into three histopathological categories: CIN1, CIN2 and CIN3. In CIN1, cells with malignant changes are limited to the superficial layer of the cervical epithelium. Most CIN1 lesions are likely to disappear without treatment. However, a small percentage may progress to high-grade CINs (i.e., CIN2 and CIN3). The risk of progression to invasive cervical cancer increases significantly with worsening CIN grades [10,11].

    Pap cytology screening for the early detection of cervical neoplasia has been successful in reducing cervical cancer incidence and mortality [12]. In unscreened populations, the risk of invasive cervical cancer occurs earlier than of most adult cancers, peaking or reaching a plateau between about 35 and 55 years of age [13]. This distribution is because cervical cancers originate mainly from HPV infections transmitted sexually in late adolescence and early adulthood [14]. HPV transmission can be reduced through the use of condoms [15]. Some studies have reported that smoking [16], multiparity [17], and long-term use of oral contraceptives [18] can double or triple the risk of precancer and cancer among women infected with carcinogenic types of HPV. There are two major kinds of anti-HPV vaccines approved for use to protect newly sexually active individuals against some of the most common HPV types and boost immunity, namely, therapeutic vaccines and prophylactic vaccines [7]. A few years after receiving a prophylactic vaccine, the individual must be revaccinated because the vaccine loses its preventive effect. Progress in the development of therapeutic vaccines for HPV has been slow [7]. In summary, there is currently no specific treatment for HPV infection [19]. There are three major treatments for cervical cancer: surgery (such as total hysterectomy and subtotal hysterectomy), radiotherapy, and chemotherapy. Among these, surgery and radiotherapy are the main treatment methods [19].

    Mathematical modeling is a useful tool for assessing the potential impact of intervention strategies against HPV spread among humans [20,21,22,23,24]. A number of authors have reported the use of mathematical modeling to evaluate the impact of HPV vaccination. Al-arydah [20] developed a two-sex, age-structured model to describe a vaccination program for the administration of an HPV vaccine. Malik et al. [11] presented an age-structured mathematical model that incorporated sex structure and Pap screening cytology. Sharomi and Malik [21] developed a two-sex HPV vaccination model to study the effect of vaccine compliance on HPV infection and cervical cancer. Omame [22] developed a two-sex deterministic model for HPV that assessed the impact of treatment and vaccination. Elbasha [23] presented a two-sex, deterministic model for assessing the potential impact of a prophylactic HPV vaccine with several properties.

    Based on the above research and understanding of HPV pathology, we develop an ordinary differential equation model with precautionary measures such as screening, which are rarely considered in previous studies, and analyze the potential effects of multiple factors on HPV transmission. The model is formulated in section 2. In section 3, the equilibria, basic reproduction number, and global stability are analyzed. We report the sensitivity analysis of the model through the partial rank correlation coefficient (PRCC) method and identify the key factors in the model in section 4. In section 5, we set the vaccination rate and screening rate as control variables and analyze an optimal control problem that minimizes vaccination and screening cost. Section 6 concludes the article. Through extensive numerical simulations with MATLAB, we obtained results to verify our conclusions.

    The total individual population at time t is divided into 10 mutually exclusive subpopulations of susceptible individuals S(t), vaccinated individuals V(t), infectious individuals without disease symptoms E(t), infectious individuals with disease symptoms H(t), individuals with persistent HPV infection P(t), CIN1 symptomatic individuals I1(t), CIN2 symptomatic individuals I2(t), CIN3 symptomatic individuals I3(t), cancer-infected individuals A(t) and recovered individuals R(t). As such, the total population is

    N(t)=V(t)+S(t)+E(t)+H(t)+P(t)+I1(t)+I2(t)+I3(t)+A(t)+R(t).

    Susceptible individuals acquire HPV infection, following effective contact with infected individuals (i.e., those in the E, H, P, I1, I2 and I3 classes) at the rate α1 as follows

    α1(t)=βαcnck(1ccca)(E+θ1H+θ2P+θ3I1+θ4I2+θ5I3)N. (1)

    It follows that the model for the transmission of HPV is given by the following system of differential equations.

    (2a)dVdt=δ2S(t)+δ4S(t)(d+δ1)V(t),(2b)dSdt=Λ+δ1V(t)(δ2+δ4+α1(t)+d)S(t),(2c)dEdt=α1(t)S(t)+δ3α1(t)R(t)+σ2H(t)(α2+cpcqγ1σ1+d)E(t),(2d)dHdt=α2E(t)+σ3P(t)(σ2+γ2σ2+α3+d)H(t),(2e)dPdt=α3H(t)+σ4I1(t)(α4+γ3σ3+σ3+d)P(t),(2f)dI1dt=α4P(t)+σ5I2(t)(α5+γ4σ4+σ4+d)I1(t),(2g)dI2dt=α5I1(t)+σ6I3(t)(σ5+γ5σ5+α6+d)I2(t),(2h)dI3dt=α6I2(t)(σ6+γ6σ6+α7+d)I3(t),(2i)dAdt=α7I3(t)(γ7+d+d1)A(t),(2j)dRdt=cpcqγ1σ1E(t)+γ2σ2H(t)+γ3σ3P(t)+γ4σ4I1(t)+γ5σ5I2(t)+γ6σ6I3(t)+γ7A(t)(δ3α1(t)+d)R(t). (2)

    Tables 1 and 2 list the associated state variables and parameters of model (2). Figure 1 shows the flow diagram of the model. We emphasize that the vaccine mentioned in model (2) is a prophylactic vaccine. In the following section, model (2) is qualitatively analyzed to derive insights into its dynamical features.

    Table 1.  Description of variables in model (2).
    Variable Description
    V(t) Vaccinated individuals
    S(t) Susceptible individuals
    E(t) Infectious individuals with no symptoms
    H(t) Infectious individuals with symptoms
    P(t) Infectious individuals with persistent infection
    I1(t) Cervical intraepithelial neoplasia grade 1 (CIN1)
    I2(t) Cervical intraepithelial neoplasia grade 2 (CIN2)
    I3(t) Cervical intraepithelial neoplasia grade 3 (CIN3)
    A(t) Cancer-infected individuals
    R(t) Recovered individuals

     | Show Table
    DownLoad: CSV
    Table 2.  Description of the parameters in model (2).
    Parameter Description
    Λ Recruitment rate into the susceptible population (per year)
    d Natural death rate (per year)
    d1 Disease-induced mortality for individuals (per year)
    α1 Effective contact rate
    δ1 Vaccine failure rate (per year)
    δ2,δ4 Vaccination rate and revaccination rate (per year)
    δ3 The modification parameter for the probability of R being infected relative to S
    cp The effect of screening by HPV testing
    cq Screening frequency (per year)
    cn Rate at which females (males) acquire new sexual partners (per year)
    ck The probability of transmitting HPV from female (male) to male (female)
    cc Condom efficacy
    ca Condom compliance (per year)
    α The negative effects of contraceptive drugs
    β The negative effects of smoking
    α2,α3,α4,α5,α6,α7 Progression rate of infectious individuals from E to H, H to P, P to I1, I1 to I2, I2 to I3, I3 to A (per year)
    σ1,σ2,σ3,σ4,σ5,σ6 Recovery rates of infectious individuals from E to R, H to E, P to H, I1 to P, I2 to I1, I3 to I2 (per year)
    γ1,γ2,γ3,γ4,γ5,γ6,γ7 Effect of drugs on infectious individuals' recovery
    θ1,θ2,θ3,θ4,θ5 Modification parameter that accounts for the infectiousness of individuals in the H, P, I1, I2, I3 classes relative to those in the E class for females (males)

     | Show Table
    DownLoad: CSV
    Figure 1.  Flow diagram of model (2).

    Model (2) is epidemiologically and mathematically well-posed in the epidemiologically valid domain

    D={(V,S,E,H,P,I1,I2,I3,A,R)R10V0,S0,E0,H0,P0,I10,I20,I30,A0,R0}.

    Theorem 3.1 Assuming that the initial condition lies in domain D, then the solutions (V,S,E,H,P,I1,I2,I3,A,R) of model (2) remain in D for all time t0. Furthermore

    lim suptN(t)Λd, with N=V+S+E+H+P+I1+I2+I3+A+R.

    Proof. We note that along the edges of D, the time derivatives all lead the solution into the invariant domain [25]

    V=0V0(2a),S=0S0(2b),E=0E0(2c),H=0H0(2d),P=0P0(2e),I1=0I10(2f),I2=0I20(2g),I3=0I30(2h),A=0A0(2i),R=0R0(2j).

    Furthermore, adding all the equations in the differential equation system of model (2) gives

    dNdt=ΛdVdSdEdHdPdI1dI2dI3dAdRd1A. (3)

    It follows from Eq (3) that

    Λ(d+d1)NdNdtΛdN.

    Therefore

    Λd+d1lim inftN(t)lim suptN(t)Λd,

    and

    lim suptN(t)Λd,

    as required.

    Model (2) is analyzed in a biologically-feasible region as follows [26]. We first show that model (2) is dissipative (i.e., all feasible solutions are uniformly bounded in a proper subset ΩR10+). Consider the region

    Ω={(V,S,E,H,P,I1,I2,I3,A,R)R10:V+S+E+H+P+I1+I2+I3+A+RΛd}.

    The following steps establish the positive invariance of Ω (i.e., solutions in Ω remain in Ωt0). It follows from Eq (3) that

    dNdtΛdN.

    A standard comparison theorem can then be used to show that

    N(t)N(0)edt+Λd(1edt).

    In particular

    N(t)Λd if N(0)Λd.

    Thus, the region Ω is positively invariant. Hence, it is sufficient to consider the dynamics of the flow generated by model (2) in Ω. In this region, the model can be considered as being epidemiologically and mathematically well-posed [27]. Thus, every solution of model (2) with initial conditions in Ω remains in Ω for all t>0. Therefore, the ω-limit sets of model (2) are contained in Ω. This result is summarized below.

    Lemma 3.1 The region Ω is positively invariant for model (2) with initial conditions in R10+.

    Model (2) has a DFE, which is obtained by setting the right-hand sides of the equations in the model to zero, given by

    ε0=(V0,S0,E0,H0,P0,I01,I02,I03,A0,R0)=(Λ(δ2+δ4)a1a2δ1(δ2+δ4),Λa1a1a2δ1(δ2+δ4),0,0,0,0,0,0,0,0). (4)

    Let X=(V,S,E,H,P,I1,I2,I3,A,R)T. Using the notation from [28], the model consists of nonnegative initial conditions together with the following system of equations:

    dXdt=Φ(X)Γ(X),

    where

    Φ(X)=(α1S000000000),Γ(X)={δ3α1Rσ2H+(α2+cpcqγ1σ1+d)Eα2Eσ3P+(σ2+γ2σ2+α3+d)Hα3Hσ4I1+(α4+γ3σ3+σ3+d)Pα4Pσ5I2+(α5+γ4σ4+σ4+d)I1α5I1σ6I3+(σ5+γ5σ5+α6+d)I2α6I2+(σ6+γ6σ6+α7+d)I3α7I3+(γ7+d+d1)AΛδ1V+(δ2+δ4+α1+d)Scpcqγ1σ1Eγ2σ2Hγ3σ3Pγ4σ4I1γ5σ5I2γ6σ6I3γ7A+(δ3α1+d)Rδ2Sδ4S+(d+δ1)V},

    and it follows that

    DΦ(X)=(F000),DΓ(X)=(VΔ0J3J4).

    The matrices F and VΔ for the new infection terms and the remaining transfer terms are respectively given by

    F=a10[1θ1θ2θ5θ4θ3000000000000000000000000000000],VΔ=[a3σ20000α2a4σ30000α3a500σ4000a8α60000σ6a7α500α40σ5a6],
    J3=[000α700cpcqγ1σ1γ2σ2γ3σ3γ6σ6γ5σ5γ4σ4a10a10θ1a10θ2a10θ5a10θ4a10θ3000000],J4=[a9000γ7d0000a2δ100δ2δ4a1],

    where

    a1=d+δ1,a2=δ2+δ4+d,a3=α2+cpcqγ1σ1+d,a4=σ2+γ2σ2+α3+d,a5=α4+γ3σ3+σ3+d,a6=α5+γ4σ4+σ4+d,a7=σ5+γ5σ5+α6+d,a8=σ6+γ6σ6+α7+d,a9=γ7+d+d1,a10=Ma1a1+δ2+δ3,M=βαcnck(1ccca).

    We obtain

    R0=ρ(FV1Δ)=a1MM1(a3D6σ2D5)(a1+δ2+δ4), (5)

    where

    M1=D6+θ1D5+θ2D4+θ3D3+θ4D2+θ5,D1=α7a9,D2=a8α6,D3=a7D2σ6α5>0,D4=a6D3σ5D2α4>0,D5=a5D4σ4D3α3>0,D6=a4D5σ3D4α2>0.

    Consequently, it follows from Theorem 2 of [28].

    Lemma 3.2 The DFE of model (2), given by (4), is locally asymptotically stable (LAS) when R0 < 1 and unstable if R0 > 1.

    The epidemiological significance of forward bifurcation is that the disease will eventually disappear if the basic reproduction number is less than one. The public health significance of backward bifurcation is that the classical requirement of R0 < 1 although necessary is no longer sufficient for effective disease control. Therefore, the presence of backward bifurcation in HPV transmission dynamics makes its effective control more difficult.

    First, the possible equilibrium solutions that model (2) can have are determined as follows. Let

    ε1=(V,S,E,H,P,I1,I2,I3,A,R),

    be any arbitrary equilibrium of model (2). Further, let

    α1=βαcnck(1ccca)(E+θ1H+θ2P+θ3I1+θ4I2+θ5I3)N, (6)

    be the associated force of infection at a steady state.

    Setting the right-hand sides of model (2) to zero (steady state) gives

    A=D1I3,I2=D2I3,I1=D3I3,P=D4I3,H=D5I3,E=D6I3,R=D7d+δ3α1I3,S=(a3D6α1σ3D7d+δ3α1σ2D5α1)I3,V=δ2+δ4a1(a3D6α1σ3D7d+δ3α1σ2D5α1)I3, (7)

    where

    D7=cpcqγ1σ1D6+γ2σ2D5+γ3σ3D4+γ4σ4D3+γ5σ5D2+γ6σ6+γ7D1.

    Substituting (7) into the expressions for α1 in (6) gives

    α1=M(D6+θ1D5+θ2D4+θ3D3+θ4D2+θ5)I3(a3D6α1σ3D7d+δ3α1σ2D5α1)(1+δ2+δ4a1)I3+D8I3+D7d+δ3α1I3, (8)

    so

    aα21+bα1+c=0, (9)

    where

    a=D8δ3a1,b=δ3(1R0)(a3D6σ2D5)(a1+δ2+δ4)+D7a1+D8a1dD7δ3(a1+δ2+δ4),c=d(1R0)(a3D6σ2D5)(a1+δ2+δ4),

    and

    D8=D6+D5+D4+D3+D2+D1+1.

    Quadratic Eq (9) can be analyzed for the possibility of multiple endemic equilibria. It is worth noting that the coefficient a is always positive, and c is positive (negative) if R0 is less than (greater than) one. Hence, the following result is established.

    Theorem 3.2 Model (2) (details in Appendix A (Table A1)) has the following.

    ⅰ. A unique endemic equilibrium if c<0R0>1;

    ⅱ. A unique endemic equilibrium if b<0, and c=0 or b24ac=0;

    ⅲ. Two endemic equilibria if c>0,b<0 and b24ac>0;

    ⅳ. No endemic equilibrium otherwise.

    Case (ⅲ) of Theorem 3.2 indicates the possibility of backward bifurcation in model (2) when R0<1. To check for this, by setting

    R1=1b2d(a3D6σ2D5)(a1+δ2+δ4),

    it can be shown that backward bifurcation occurs for values of R1<R0<1. This phenomenon is illustrated by simulating model (2). The parameter values are presented in Table 3. Let M[0.35,0.5]. It should be mentioned that the aforementioned parameter values may not all be epidemiologically realistic.

    Table 3.  Parameter values used in Figure 2 (A: Assumed).
    Parameter Value Source Parameter Value Source Parameter Value Source
    Λ 288802 [22] α3 0.005 [11] γ1 1.5 A
    d 0.0162 [22] α4 0.1 [11] γ2 1.5 A
    d1 0.01 [11] α5 0.02 [11] γ3 1.2 A
    δ1 0.1 [11] α6 0.04 [11] γ4 1.1 A
    δ2 0.87 [22] α7 0.08 [11] γ5 1.05 A
    δ3 0.3 [22] σ1 0.99 [11] γ6 1.03 A
    δ4 0.27 A σ2 9e-4 [22] γ7 1.01 A
    cp 0.9 A σ3 0.5 [22] θ4 0.6 A
    cq 0.4 A σ4 1.9e-7 A θ5 0.5 A
    σ5 1.9e-7 A θ1 1 A θ2 0.8 [22]
    α2 0.5 [22] σ6 1.9e-7 A θ3 0.7 A

     | Show Table
    DownLoad: CSV
    Figure 2.  Backward bifurcation diagram of model (2).

    The associated backward bifurcation diagram, depicted in Figure 2, shows that the model has a DFE (corresponding to Figure 3) and two endemic equilibria: One of the endemic equilibria is LAS (corresponding to Figure 4a); the other is unstable (a saddle); and the disease-free equilibrium is LAS. This clearly shows the coexistence of two stable equilibria when R0<1, confirming that the model exhibits backward bifurcation for R1<R0<1. This result is summarized below for model (2) (a more rigorous proof of the backward bifurcation phenomenon of the model, using the center manifold theory is given in Appendix B).

    Figure 3.  Variation in population with R0=0.6659 and R1=0.9976.
    Figure 4.  Variation in population with (a) R0=0.9988, R1=0.8144; (b) R0=2.2196.

    Theorem 3.3 Model (2) exhibits backward bifurcation when Case (ⅲ) of Theorem 3.2 holds and R1<R0<1.

    Consider model (2) with perfect protection after recovery (that is, δ3=0). In such a case, the basic reproduction number is R0=R0|δ3=0. It follows from Eq (9) that if δ3=0, the coefficients a=0 and b>0, so quadratic Eq (9) reduces to a linear equation in α1 (with α1=c/b). In this case, model (2) has a unique endemic equilibrium if c<0 (i.e., R0>1), ruling out backward bifurcation in the model for this case (the presence of two endemic equilibria when R0<1 is necessary for the existence of backward bifurcation). Furthermore, it follows that c=0 when R0=1. Thus, in such a case (with a=c=0), quadratic Eq (9) has only the trivial solution α1=0 (which corresponds to the DFE ε0). This result is summarized below.

    Lemma 3.3 Consider the case where the protection after recovery is perfect (δ3=0). Model (2) has a unique endemic equilibrium whenever R0>1 and no endemic equilibrium otherwise.

    Theorem 3.4 In the first quadrant, there is no limit cycle in model (2).

    Proof We consider the Dulac function as B(S,E)=1SE. Let

    Q=E+θ1H+θ2P+θ3I1+θ4I2+θ5I3.

    Hence QE and NR. Therefore,

    MQNMQRN20,EME2NMQE2N0.

    Then,

    Θ=(BV)V+(BS)S+(BE)E+(BH)H+(BP)P+(BI1)I1+(BI2)I2+(BI3)I3+(BA)A+(BR)R=a1SEΛ+δ1VS2E+MQEN2+[EME2NMQEN2MQE2N]+δ3RS[EME2NMQEN2MQE2N]σ2HSE2a4SEa5SEa5SEa6SEa7SEa8SEa9SEdSEδ3SE(MQNMQRN2)<0.

    Therefore, by the Dulac−Bendixson theorem [29], there is no periodic orbit for model (2). Moreover, ε0 is the unique positive equilibrium point in R10+ if δ3=0, and it is also locally asymptotically stable for R0<1. Hence, every positive solution actually approaches ε0. Thus, ε0 is globally asymptotically stable if δ3=0 and R0<1.

    In this section, we performed a numerical simulation to enhance the understanding of model (2).

    To examine the possible impact of interventions on disease infections we plot the number of infected individuals (E) with various vaccination rates and revaccination rates.

    This analysis shows that an increasing vaccination rate persistently decreases the peak value, as shown in Figure 5. Increasing the vaccination rate δ2 by 1.75 times (increase from 0.4 to 0.7) or 1.43 times (increased from 0.7 to 1) will lead to a reduction in the peak value in the number of E by 20.21% or by 15.67%, respectively. In addition, the peak value of the number of people infected with δ2=1 decreased by 43.82% compared with the number of people infected with δ2=0.

    Figure 5.  Variation in population E with different parameters (a) δ2; (b) δ2+δ4.

    On the premise that the vaccine's protective effect will end after a few years, we consider the situation of vaccination and revaccination. Figure 5b indicates that increasing δ2 and δ4 from 0 to 0.4 will lead to a reduction in the peak value in the number of E by 34.16%. In addition, the peak value of the number of people infected with δ2=δ4=0.7 decreased by 100% compared with the number of people infected with δ2=δ4=0.

    To identify the factors associated with a certain intervention that markedly affect the rate of new infections, we performed sensitivity analysis of the basic reproduction number.

    LHS belongs to the MC class of sampling methods; it was introduced by Mckay et al. [30]. LHS allows an unbiased estimate of the average model output and has the advantage that it requires fewer samples than simple random sampling to achieve the same accuracy. For nonlinear but monotonic relationships between outputs and inputs, measures that work well are based on rank transforms such as the partial rank correlation coefficient, and standardized rank regression coefficient.

    Model (2) has 39 parameters. To identify the key factors, following [30], we performed a Latin hypercube sampling on the parameters that appear in R0 and calculated the PRCC. The parameters of the model were set as input variables, and R0 was the output variable. Generally, in PRCC analysis, the parameters with large PRCC values and corresponding small p values are deemed to be the most influential parameters in the model.

    Detailed inspection of Table C1 (Appendix C) and Figure 6 indicates that in terms of reducing the value of R0, except σ3 (control the disease and reduce the number of persistent infections) and d, the vaccination rate δ2 is the most sensitive parameter with a leading PRCC value, followed by γ2,γ3,σ2,δ4. This implies that enhancing the vaccination rate is the most effective intervention for lowering HPV new infections. Moreover, in the treatment of patients in stages H, P, I1, I2 and I3, the effect of treatments γ2,γ3,γ4,γ5 and γ6 on R0 decreases successively. That is, the same treatment intervention is more effective in the earlier stages. This means that more attention should be paid to patients in the early stages of infection. As asymptomatic patients are unable to diagnose themselves, regular screening for HPV should be strengthened. Smoking, overuse of contraceptive drugs, and unsafe sexual life will increase the value of R0, thus promoting the spread of HPV.

    Figure 6.  Significance test of model parameters and PRCC results for R0.

    In this section, an optimal control model for the transmission dynamics of HPV is formulated by extending model (2) to include control functions. Our goal here is to study the optimal control strategies to curtail the epidemic and minimize cost.

    The optimal vaccination and screening strategy can be formulated as the following optimal control problem (P) with inequality constraints and free terminal states defined over the prescribed interval [0,tf] [31]:

    minJ=tf0(C1E2+C2H2+C3P2+C4I21+C5I22+C6I23+B1u21+B2u22)dt,s.t.V=u1(t)S+δ4S(d+δ1)V,S=Λ+δ1V(u1(t)+δ4+α1+d)S,E=α1S+δ3α1R+σ2H(α2+cpu2(t)γ1σ1+d)E,H=α2E+σ3P(σ2+γ2σ2+α3+d)H,P=α3H+σ4I1(α4+γ3σ3+σ3+d)P,I1=α4P+σ5I2(α5+γ4σ4+σ4+d)I1,I2=α5I1+σ6I3(σ5+γ5σ5+α6+d)I2,I3=α6I2(σ6+γ6σ6+α7+d)I3,A=α7I3(γ7+d+d1)A,R=cpu2(t)γ1σ1E+γ2σ2H+γ3σ3P+γ4σ4I1+γ5σ5I2+γ6σ6I3+γ7A(δ3α1+d)R,V(0)=Vs,S(0)=Ss,E(0)=Es,H(0)=Hs,P(0)=Ps,I1(0)=I1s,I2(0)=I2s,I3(0)=I3s,A(0)=As,R(0)=Rs,0u1(t)u1max,0u2(t)u2max, (10)

    where tfR+ is the fixed terminal time, the coefficients C1,C2,C3,C4,C5,C6,B1 and B2 represent the corresponding weight constants, and these weights are balancing cost factors related to the size and importance of the parts of the objective function. The control function u1(t) is the fraction of the population of susceptible individuals who enters the vaccination compartment. The control function u2(t) is the fraction of the population of infectious individuals with no symptoms who undergo HPV screening, and they are Lebesgue integrable.

    The inequality constraints in problem (P) can be transformed into equality ones with the help of some non-negative parametric parameters, that is, ηi(i=1,2,3,4), as

    {u1+η1=0,u1u1max+η2=0,u2+η3=0,u2u2max+η4=0. (11)

    Hence, the Hamiltonian function for problem (P) is obtained as follows:

    HΔ=C1E2+C2H2+C3P2+C4I21+C5I22+C6I23+B1u21+B2u22+λV[u1S+δ4S(d+δ1)V]+λS[Λ+δ1V(u1+δ4+α1+d)S]+λE[α1S+δ3α1R+σ2H(α2+cpu2γ1σ1+d)E]+λH[α2E+σ3P(σ2+γ2σ2+α3+d)H]+λI1[α4P+σ5I2(α5+γ4σ4+σ4+d)I1]+λI2[α5I1+σ6I3(σ5+γ5σ5+α6+d)I2]+λI3[α6I2(σ6+γ6σ6+α7+d)I3]+λA[α7I3(γ7+d+d1)A]+λR[cpu2γ1σ1E+γ2σ2H+γ3σ3P+γ4σ4I1+γ5σ5I2+γ6σ6I3+γ7A(δ3α1+d)R]+μ1(u1+η1)+μ2(u1u1max+η2)+μ3(u2+η3)+μ4(u2u2max+η4), (12)

    where λ=[λV,λS,λE,λH,λP,λI1,λI2,λI3,λA,λR]T are adjoint variables, and μ=[μ1,μ2,μ3,μ4]T are non-negative penalty multipliers [32].

    Theorem 5.1 There exists an optimal control (u1(t),u2(t)) and corresponding solution V, S, E, H, P, I1,I2,I3,A, and R that minimize J(u1(t),u2(t)) over Ω. Furthermore, there exist adjoint functions λV,λS,λE,λH,λP,λI1,λI2,λI3,λA and λR, such that

    λV=λS(δ1+α1SN)+λEα1S+α1δ3RN+λVa1λRα1δ3RN,
    λS=λS(d+α1+δ4+u1α1SN)+λE(α1S+α1δ3RNα1)λRα1δ3RNλV(u1+δ4),
    λE=2C1EλHα2+λE(a3+α1S+δ3RMSMRδ3N)+λSS(Mα1)N
    λR(cpu2γ1σ1+α1δ3Rα1δ3N),
    λH=2C2HλPα3λE(σ2+MSθ1+MRδ3θ1α1δ3Rα1SN)+λHa4+λSS(Mθ1α1)N
    λR(γ2σ2+α1δ3RMRδ3θ1N),
    λP=2C3PλI1α4λHσ3+λPa5λE(σ2+MSθ2+MRδ3θ2α1δ3Rα1SN)
    +λSS(Mθ2α1)NλR(γ3σ3+α1δ3RMRδ3θ2N),
    λI1=2C4I1λI2α5λPσ4+λSS(Mθ3α1)NλEMSθ3+MRδ3θ3α1δ3Rα1SN
    +λI1a6λR(γ4σ4+α1δ3RMRδ3θ3N),
    λI2=2C5I2λI3α6λI1σ5+λI2a7λEMSθ4+MRδ3θ4α1δ3Rα1SN
    λR(γ5σ5+α1δ3RMRδ3θ4N)+λSS(Mθ4α1)N,
    λI3=2C6I3λAα7λI2σ6+λI3a8λEMSθ5+MRδ3θ5α1δ3Rα1SN
    λR(γ6σ6+α1δ3RMRδ3θ5N)+λSS(Mθ5α1)N,
    λA=λR(γ7+α1δ3RN)+λEα1δ3R+α1SN+λAa9λSSα1N,
    λR=λR(d+δ3α1α1δ3RN)λSSα1N+λE(α1δ3R+α1SNδ3α1), (13)

    with transversality conditions

    λV(tf)=λS(tf)=L=λA(tf)=λR(tf)=0. (14)

    The following characterization holds

    {u1(t)=max{0,min{S(λSλV)2B1,u1max}},u2(t)=max{0,min{Ecpγ1σ1(λEλR)2B2,u2max}}. (15)

    Proof. The existence of an optimal control can be obtained owing to the convexity of the integrand of J(u1(t),u2(t)) with respect to (u1(t),u2(t)) [33], a priori boundedness of the state solutions, and the Lipschitz property of the state system with respect to the state variables.

    By Pontryagin's maximum principle [34], the optimal conditions with respect to the state, costate, and parametric variables result in a two-point boundary value problem coupled with a nonlinear complementarity problem as follows:

    λV=HΔV,λS=HΔS,L,λA=HΔA,λR=HΔR, (16)

    and

    λV(tf)=λS(tf)=L=λA(tf)=λR(tf)=0,

    evaluated at the optimal control and corresponding states results in the stated adjoint system (13) with transversality (14).

    The optimality conditions with respect to the control variables are

    HΔu1=0,HΔu2=0. (17)

    By solving Eq (17), the optimal control can be expressed as

    u1(t)=S(λSλV)+μ1μ22B1.

    To determine an explicit expression of the optimal control without μ1, we consider the following three cases:

    ⅰ. On the set {t0<u1<u1max}, we have μ1=μ2=0. Hence, u1(t)=S(λSλV)2B1.

    ⅱ. On the set {tu1=u1max}, we have μ2=0. Hence, u1(t)=u1max=S(λSλV)+μ12B1. As μ10, it is determined that u1maxS(λSλV)2B1.

    ⅲ. On the set {tu1=0}, we have μ1=0. Hence, u1(t)=0=S(λSλV)μ22B1. As μ20, it is determined that S(λSλV)2B10.

    Combining the above three cases, the optimal control u1 is characterized as

    u1(t)=max{0,min{S(λSλV)2B1,u1max}}. (18)

    Using similar arguments, we can characterize the optimal control u2 as

    u2(t)=max{0,min{Ecpγ1σ1(λEλR)2B2,u2max}}. (19)

    An analytical expression of the optimal vaccination rate and screening rate was derived in Eq (15). However, an effective algorithm is still required to solve the nonlinear constrained optimal control problem numerically. Based on the generating function method, Peng et al. developed a series of symplectic methods for nonlinear optimal control problems [35,36,37,38]. Such symplectic methods have good precision and efficiency because of the structure-preserved property. Recently, Wang et al. improved the symplectic methods by incorporating the local pseudospectral discretization scheme [39,40,41,42]. Such symplectic pseudospectral methods (SPMs) have been successfully applied to solve optimal control problems in various mechanical systems [43,44]. In this paper, the SPM developed in [45] was adopted.

    In the following simulation, the weights in the objective function (meaning the minimization of the number of patients at each stage has different importance) are {C_1} = 4.5e - 2, \quad {C_2} = 1e - 7, \; {C_3} = 1e - 4, \quad {C_4} = 1e - 5, \quad {C_5} = 2e - 4, and C_{6} = 1 e-4. Let M = 1. The initial values for the states and other parameters are listed in Table 4. Unless otherwise stated, the parameters used in each case were as listed in Table 3.

    Table 4.  Parameter values used in Figure 7.
    Parameter Value Parameter Value
    Vs 3.2607e6 I2s 4.15e4
    Ss 3.2212e5 I3s 1.68e4
    Es 4.2204e4 As 1.29e3
    Hs 1.1162e7 Rs 6.3e3
    Ps 4e4 tf 50
    I1s 1.15e5 u_{1 \max } 1
    u_{2 \max } 2

     | Show Table
    DownLoad: CSV
    Figure 7.  Simulations of model (2). Dashed lines: Populations with optimal control. Solid lines: Populations without control. Parameter values are B_{1} = 8.3 e 8 and B_{2} = 4 e 8.

    The controlled solutions together with the solutions for the uncontrolled case are presented in Figure 7. It can be seen that the control strategy is effective. Vaccinated individuals increase steadily and reach almost 400% at the terminal end. Susceptible individuals keep increasing and then stabilize during the whole period. The number of infected individuals decreases significantly when optimal control strategies are used compared to the number in the absence of control strategies.

    We considered another set of weights, the simulation results are shown in Figure 8. A higher focus on the control strategies leads to a drop in the importance of the vaccination and screening strategies. As the number of asymptomatic individuals depends on the immunity of the susceptible individuals and the protection of the susceptible population, we should consider strengthening their immunity or implement regular cost-effective screening to control HPV transmission.

    Figure 8.  The different strategies of u_{1}(t) and u_{2}(t) are plotted for B_{1} = 4 e 9 and B_{2} = 3 e 9. Other parameter values are the same as those in Figure 7.

    The human papillomavirus is among the most common sexually transmitted infections. Following infection, cervical carcinogenesis is a complex stepwise process characterized by slow progression. According to the known pathology, we represented the CIN stages with three corresponding components in the model. Our model accounted for the fact that preventive vaccines become ineffective over time. We derived three types of equilibria and their conditions of existence, analyzed the stability of the equilibria, and characterized the threshold condition as backward bifurcation for the stable fixed points. We also obtained the conditions for the elimination of the disease. We found that the possibility of HPV transmission to lead to endemic disease can be reduced by strengthening the protection after cure. We then simulated and compared practical mitigation strategies and performed sensitivity analysis to illustrate the key factors for the threshold condition. The results show that increasing the vaccination rate is the most effective way to reduce the basic reproduction number. The effect of optimal control was illustrated numerically, and a comparison of HPV infection was presented under different control strategies.

    This work was supported by the Fundamental Research Funds for the Central Universities (31920200037; 31920200070), the Research Fund for Humanities and Social Sciences of the Ministry of Education(20XJAZH006), the Program for Young Talent of State Ethnic Affairs Commission of China (No. [2014]121), the Innovation Team of Intelligent Computing and Dynamical System Analysis and Application.

    The authors declare that no conflict of interest.

    Table A1.  A detailed explanation of Theorem 3.2.
    a c b Results
    a>0 c>0 b>0 Two negative points
    b=0 No equilibrium points
    b < 0 Two endemic equilibria if b^{2}-4 a c>0
    c=0 b>0 A negative point and a DFE
    b=0 Two DFE
    b < 0 A DFE and an EEP (or b^{2}-4 a c=0)
    c < 0 b>0 A negative point and an EEP
    b=0 A negative point and an EEP
    b < 0 A negative point and an EEP

     | Show Table
    DownLoad: CSV

    Here, we explore the existence of backward bifurcation using the center manifold theory [46,47]. To apply this theory, it is necessary to carry out the following change of variables.

    \text { Let } E = x_{1}, \quad H = x_{2}, \quad P = x_{3}, \quad I_{3} = x_{4}, \quad I_{2} = x_{5}, \quad I_{1} = x_{6}, \quad A = x_{7}, \quad R = x_{8}, \quad S = x_{9},

    V = x_{10}, so that

    N = \sum\limits_{i = 1}^{10} x_{i}.

    Further, using the vector notation

    X_{1} = \left(x_{1}, x_{2}, x_{3}, x_{4}, x_{5}, x_{6}, x_{7}, x_{8}, x_{9}, x_{10}\right)^{T}.

    Model (2) can be rewritten in the form

    \frac{d X_{1}}{d t} = F = \left(f_{1}, f_{2}, f_{3}, f_{4}, f_{5}, f_{6}, f_{7}, f_{8}, f_{9}, f_{10}\right)^{T},

    as follows:

    \left\{ \begin{array}{l} \frac{{d{x_1}}}{{dt}} = {\alpha _1}(t){x_9}(t) + {\delta _3}{\alpha _1}(t){x_7}(t) + {\sigma _2}{x_3}(t) - \left( {{\alpha _2} + {c_p}{c_q}{\gamma _1}{\sigma _1} + d} \right){x_1}(t), \\ \frac{{d{x_2}}}{{dt}} = {\alpha _2}{x_1}(t) + {\sigma _3}{x_3}(t) - \left( {{\sigma _2} + {\gamma _2}{\sigma _2} + {\alpha _3} + d} \right){x_2}(t), \\ \frac{{d{x_3}}}{{dt}} = {\alpha _3}{x_2}(t) + {\sigma _4}{x_6}(t) - \left( {{\alpha _4} + {\gamma _3}{\sigma _3} + {\sigma _3} + d} \right){x_3}(t), \\ \frac{{d{x_4}}}{{dt}} = {\alpha _6}{x_5}(t) - \left( {{\sigma _6} + {\gamma _6}{\sigma _6} + {\alpha _7} + d} \right){x_4}(t), \\ \frac{{d{x_5}}}{{dt}} = {\alpha _5}{x_6} + {\sigma _6}{x_4}(t) - \left( {{\sigma _5} + {\gamma _5}{\sigma _5} + {\alpha _6} + d} \right){x_5}(t), \\ \frac{{d{x_6}}}{{dt}} = {\alpha _4}{x_3}(t) + {\sigma _5}{x_5}(t) - \left( {{\alpha _5} + {\gamma _4}{\sigma _4} + {\sigma _4} + d} \right){x_6}(t), \\ \frac{{d{x_7}}}{{dt}} = {\alpha _7}{x_4}(t) - \left( {{\gamma _7} + d + {d_1}} \right){x_7}(t), \\ \frac{{d{x_8}}}{{dt}} = {c_p}{c_q}{\gamma _1}{\sigma _1}{x_1}(t) + {\gamma _2}{\sigma _2}{x_2}(t) + {\gamma _3}{\sigma _3}{x_3}(t) + {\gamma _4}{\sigma _4}{x_5}(t) + {\gamma _5}{\sigma _5}{x_5}(t) + {\gamma _6}{\sigma _6}{x_4}(t) + {\gamma _7}{x_7}(t)\\ \;\;\;\;\;\;\; - \left( {{\delta _3}{\alpha _1}(t) + d} \right){x_8}(t), \\ \frac{{d{x_9}}}{{dt}} = \Lambda + {\delta _1}{x_{10}}(t) - \left( {{\delta _2} + {\delta _4} + {\alpha _1}(t) + d} \right){x_9}(t), \\ \frac{{d{x_{10}}}}{{dt}} = {\delta _2}{x_9}(t) + {\delta _4}{x_9}(t) - \left( {d + {\delta _1}} \right){x_{10}}(t), \end{array} \right. (B.1)

    with

    \alpha_{1}(t) = \frac{\beta \alpha c_{n} c_{k}\left(1-c_{c} c_{a}\right)\left(x_{1}+\theta_{1} x_{2}+\theta_{2} x_{3}+\theta_{3} x_{6}+\theta_{4} x_{5}+\theta_{5} x_{4}\right)}{N}.

    Consider the case when R_{0} = 1. Suppose, further, that {\delta _2} is chosen as a bifurcation parameter.

    Solving for \delta_{2} = \delta_{2^{*}} from {R_0} gives

    \delta_{2^{*}} = \frac{a_{1} M M_{1}}{\left(a_{3} D_{6}-\sigma_{2} D_{5}\right)}-a_{1}-\delta_{4}.

    The Jacobian of model (B.1) evaluated at the DFE is given as

    J\left(\varepsilon_{0}\right) = \left[\begin{array}{cc} J_{11} & J_{12} \\ J_{21} & J_{22} \end{array}\right],

    where

    {J_{11}} = \left[ {\begin{array}{*{20}{c}} {{a_{10}} - {a_3}}&{{a_{10}}{\theta _1} + {\sigma _2}}&{{a_{10}}{\theta _2}}&{{a_{10}}{\theta _5}}&{{a_{10}}{\theta _4}}&{{a_{10}}{\theta _3}}\\ {{\alpha _2}}&{ - {a_4}}&{{\sigma _3}}&0&0&0\\ 0&{{\alpha _3}}&{ - {a_5}}&0&0&{{\sigma _4}}\\ 0&0&0&{ - {a_8}}&{{\alpha _6}}&0\\ 0&0&0&{{\sigma _6}}&{ - {a_7}}&{{\alpha _5}}\\ 0&0&{{\alpha _4}}&0&{{\sigma _5}}&{ - {a_6}} \end{array}} \right], \;{J_{12}} = {[0]_{6 \times 4}},
    {J_{21}} = \left[ {\begin{array}{*{20}{c}} 0&0&0&{{\alpha _7}}&0&0\\ {{c_p}{c_q}{\gamma _1}{\sigma _1}}&{{\gamma _2}{\sigma _2}}&{{\gamma _3}{\sigma _3}}&{{\gamma _6}{\sigma _6}}&{{\gamma _5}{\sigma _5}}&{{\gamma _4}{\sigma _4}}\\ { - {a_{10}}}&0&{ - {a_{10}}{\theta _1}}&{ - {a_{10}}{\theta _4}}&{ - {a_{10}}{\theta _3}}&{ - {a_{10}}{\theta _2}}\\ 0&0&0&0&0&0 \end{array}} \right],
    {J_{22}} = \left[ {\begin{array}{*{20}{c}} { - {a_9}}&0&0&0\\ {{\gamma _7}}&{ - d}&0&0\\ 0&0&{ - {a_2}}&{{\delta _1}}\\ 0&0&{{\delta _2} + {\delta _4}}&{ - {a_1}} \end{array}} \right].

    It is easy to verify that the transformed model (B.1), with \delta_{2} = \delta_{2^{*}}, has a hyperbolic equilibrium point (i.e., the linearized system has a simple eigenvalue with zero real part, and all other eigenvalues have negative real parts). Hence, the center manifold theory can be used to analyse the dynamics of model (B.1) near \delta_{2} = \delta_{2^{*}}.

    It can be shown that the Jacobian of model (B.1) at \delta_{2} = \delta_{2^{*}} has a right eigenvector (associated with the zero eigenvalue) given by w = \left(w_{1}, w_{2}, w_{3}, w_{4}, w_{5}, w_{6}, w_{7}, w_{8}, w_{9}, w_{10}\right)^{T}, where

    w_{1} = \frac{a_{4} D_{5}-\sigma_{3} D_{4}}{\alpha_{2}} w_{4} = D_{6} w_{4} \gt 0, \quad w_{2} = \frac{a_{5} D_{4}-\sigma_{4} D_{3}}{\alpha_{3}} w_{4} = D_{5} w_{4} \gt 0,
    w_{3} = \frac{a_{6} D_{3}-\sigma_{5} D_{2}}{\alpha_{4}} w_{4} = D_{4} w_{4} \gt 0, \quad w_{4} = w_{4} \gt 0, \quad w_{5} = \frac{a_{8}}{\alpha_{6}} w_{4} = D_{2} w_{4} \gt 0,
    w_{6} = \frac{a_{7} D_{2}-\sigma_{6}}{\alpha_{5}} w_{4} = D_{3} w_{4} \gt 0, \quad w_{7} = \frac{\alpha_{7}}{a_{9}} w_{4} = D_{1} w_{4} \gt 0,
    {w_8} = \frac{{{c_p}{c_q}{\gamma _1}{\sigma _1}{D_6} + {\gamma _2}{\sigma _2}{D_5} + {\gamma _3}{\sigma _3}{D_4} + {\gamma _6}{\sigma _6} + {\gamma _5}{\sigma _5}{D_2} + {\gamma _4}{\sigma _4}{D_3} + {\gamma _7}{D_1}}}{d}{w_4} \gt 0,
    w_{9} = \frac{a_{1}}{\delta_{2}+\delta_{4}} w_{10} = G_{1} G_{2} w_{4} = G_{3} w_{4}, \quad w_{10} = \frac{\sigma_{2} D_{5}-a_{3} D_{6}}{a_{2} G_{1}-\delta_{1}} w_{4} = G_{2} w_{4}.

    The components of the left eigenvector of \left.J_{\varepsilon_{0}}\right|_{\delta_{2} = \delta_{2^{*}}}, \quad v = \left(v_{1}, v_{2}, v_{3}, v_{4}, v_{5}, v_{6}, v_{7}, v_{8}, v_{9}, v_{10}\right), satisfying v \cdot w = 1 are

    v_{1} \gt 0, \quad v_{2} = \frac{a_{3}}{\alpha_{2}} v_{1}, \quad v_{3} = \frac{a_{3} a_{4}-\alpha_{2} \sigma_{2}}{\alpha_{2} \alpha_{3}} v_{1} \gt 0, \quad v_{4} = \frac{\sigma_{6}}{a_{8}} v_{5},
    v_{5} = \frac{a_{3} a_{4} a_{5} a_{6}-a_{5} a_{6} \sigma_{2} \alpha_{2}-a_{3} a_{6} \sigma_{3} \alpha_{3}-a_{3} a_{4} \sigma_{4} \alpha_{4}+\sigma_{4} \alpha_{4} \sigma_{2} \alpha_{2}}{\alpha_{2} \alpha_{3} \alpha_{4} \alpha_{5}} \gt 0,
    v_{6} = \frac{a_{3} a_{4} a_{5}-a_{5} \sigma_{2} \alpha_{2}-a_{3} \sigma_{3} \alpha_{3}}{\alpha_{2} \alpha_{3} \alpha_{4}} \gt 0, \quad v_{7} = v_{8} = 0, \quad v_{9} = v_{1}, \quad v_{10} = \frac{\delta_{1}}{a_{1}} v_{1}.

    It follows from [26]:

    a = \sum\limits_{k, i, j = 1}^{10} v_{k} \omega_{i} \omega_{j} \frac{\partial^{2} f_{k}}{\partial x_{i} \partial x_{j}}(0, 0), \quad b = \sum\limits_{k, i = 1}^{10} v_{k} \omega_{i} \frac{\partial^{2} f_{k}}{\partial x_{i} \partial \delta_{2^{*}}}(0, 0),

    are computed to be

    \begin{array}{l} a = \sum\limits_{k, i, j = 1}^{10} {{v_k}} {w_i}{w_j}\frac{{{\partial ^2}{f_k}}}{{\partial {x_i}\partial {x_j}}}(0, 0)\\ \;\; = \frac{{2{v_1}{w_8}M{\delta _3}}}{{{S^0} + {V^0}}}\left( {{w_1} + {\theta _1}{w_2} + {\theta _2}{w_3} + {\theta _5}{w_4} + {\theta _3}{w_6}} \right)\\ \;\;{\rm{ \gt 0}}, \end{array} (B.2)
    \begin{array}{l} b = \sum\limits_{k, i = 1}^{10} {{v_k}} {w_i}\frac{{{\partial ^2}{f_k}}}{{\partial {x_i}\partial {\delta _{{2^*}}}}}(0, 0)\\ = {v_9}{w_9}\frac{{{\partial ^2}{f_9}}}{{\partial {x_9}\partial {\delta _{{2^*}}}}}(0, 0) + {v_{10}}{w_9}\frac{{{\partial ^2}{f_{10}}}}{{\partial {x_9}\partial {\delta _{{2^*}}}}}(0, 0)\\ = - {v_1}\frac{{{\sigma _2}{D_5} - {a_3}{D_6}}}{{{a_2}{G_1} - {\delta _1}}} \times \frac{{{a_1}}}{{{\delta _2} + {\delta _4}}}{w_4} + \frac{{{\delta _1}}}{{{a_1}}}{v_1}\frac{{{\sigma _2}{D_5} - {a_3}{D_6}}}{{{a_2}{G_1} - {\delta _1}}} \times \frac{{{a_1}}}{{{\delta _2} + {\delta _4}}}{w_4}\\ = \left( {\frac{{{\delta _1}}}{{{a_1}}} - 1} \right){v_1}\frac{{{\sigma _2}{D_5} - {a_3}{D_6}}}{{{a_2}\left( {\frac{{{a_1}}}{{{\delta _2} + {\delta _4}}}} \right) - {\delta _1}}} \times \frac{{{a_1}}}{{{\delta _2} + {\delta _4}}}{w_4}\\ {\rm{ \gt 0}}. \end{array} (B.3)

    Thus, we have made the following conclusions

    Theorem A.1 Model (B.1) (or, equivalently, model (2)) undergoes a backward bifurcation at R_{0} = 1 if all parameters are positive.

    Table C1.  PRCC values of R0 with corresponding values of p (significant for p ≤ 0.01).
    Parameter PRCC p values
    θ1 0.1079 1.3332e-06
    θ2 0.0619 0.0056
    θ3 0.0375 0.0939
    θ4 0.0363 0.1050
    θ5 0.0130 0.5625
    δ1 0.2219 1.0079e-23
    δ2 −0.0972 1.3420e-05
    δ3 0.0139 0.5334
    δ4 −0.0503 0.0246
    σ1 −0.0474 0.0342
    σ2 −0.0907 4.8958e-05
    σ3 −0.1749 3.3636e-15
    σ4 −0.0631 0.0047
    σ5 −0.0038 0.8662
    σ6 0.0162 0.4703
    γ1 −0.0491 0.0280
    γ2 −0.1101 7.9369e-07
    γ3 −0.1100 8.1799e-07
    γ4 −0.0857 0.0001
    γ5 −0.0117 0.6015
    γ6 −0.0179 0.4227
    γ7 −0.0254 0.2562
    α2 0.1070 1.6317e-06
    α3 0.0391 0.0806
    α4 0.0389 0.0823
    α5 0.0347 0.1209
    α6 −0.0100 0.6578
    α7 −0.0331 0.1389
    \Lambda 0.0106 0.6363
    d1 0.0018 0.9352
    d −0.2728 1.7789e-35
    cp −0.1683 3.5473e-14
    cq −0.1044 2.9039e-06
    cc −0.0685 0.0022
    ca −0.0648 0.0037
    cn 0.3778 7.0664e-69
    ck 0.2383 3.1497e-27
    \alpha 0.2481 1.9591e-29
    \beta 0.6344 1.0688e-225

     | Show Table
    DownLoad: CSV


    [1] H. Gao, K. Mangion, D. Carrick, D. Husmeier, X. Y. Luo, C. Berry, Estimating prognosis in patients with acute myocardial infarction using personalized computational heart models, Sci. Rep., 7 (2017), 13527. https://doi.org/10.5465/AMBPP.2017.13527abstract doi: 10.5465/AMBPP.2017.13527abstract
    [2] K. Mangion, H. Gao, D. Husmeier, X. Y. Luo, C. Berry, Advances in computational modelling for personalised medicine after myocardial infarction, Heart (British Cardiac Society), 104 (2018), 550–557. https://doi.org/10.1136/heartjnl-2017-311449 doi: 10.1136/heartjnl-2017-311449
    [3] D. D. Streeter, W. T. Hanna, Engineering mechanics for successive states in canine left ventricular myocardium, Circ. Res., 33 (1973), 639–55. https://doi.org/10.1161/01.RES.33.6.639 doi: 10.1161/01.RES.33.6.639
    [4] S. F. Pravdin, V. I. Berdyshev, A. V. Panfilov, L. B. Katsnelson, V. S. Markhasin, Mathematical model of the anatomy and fibre orientation field of the left ventricle of the heart, Biomed. Eng. Online, 12 (2013), 54. https://doi.org/10.1186/1475-925X-12-54 doi: 10.1186/1475-925X-12-54
    [5] P. D. Achille, A. Harouni, S. Khamzin, O. Solovyova, J. J. Rice, Gaussian process regressions for inverse problems and parameter searches in models of ventricular mechanics, Front. Physiol., 9 (2018), 1002. https://doi.org/10.3389/fphys.2018.01002 doi: 10.3389/fphys.2018.01002
    [6] P. Gould, D. Ghista, L. Brombolich, I. Mirsky, In vivo stresses in the human left ventricular wall: Analysis accounting for the irregular 3-dimensional geometry and comparison with idealised geometry analyses, J. Biomech., 5 (1972), 521–539. https://doi.org/10.1016/0021-9290(72)90009-7 doi: 10.1016/0021-9290(72)90009-7
    [7] M. Sermesant, P. Moireau, O. Camara, J. S. Marie, R. Andriantsimiavona, R. Cimrman, et al., Cardiac function estimation from MRI using a heart model and data assimilation: advances and difficulties, Med. Image Anal., 10 (2006), 642–656. https://doi.org/10.1016/j.media.2006.04.002 doi: 10.1016/j.media.2006.04.002
    [8] C. S. Peskin, Flow patterns around heart valves: A numerical method, J. Comput. Phys., 10 (1972), 252–271. https://doi.org/10.1016/0021-9991(72)90065-4 doi: 10.1016/0021-9991(72)90065-4
    [9] L. Cai, H. Gao, X. Y. Luo, Y. F. Nie, Multi-scale modelling of the human left ventricle, Sci. Sin., 45 (2015), 024702. https://doi.org/10.1360/SSPMA2013-00100 doi: 10.1360/SSPMA2013-00100
    [10] B. E. Griffith, X. Y. Luo, Hybrid finite difference/finite element immersed boundary method, Int. J. Numer. Methods Biomed. Eng., 33 (2017), 101002. https://doi.org/10.1002/cnm.2888 doi: 10.1002/cnm.2888
    [11] H. Gao, H. M. Wang, C. Berry, X. Y. Luo, B. E. Griffith, Quasi-static imaged-based immersed boundary-finite element model of human left ventricle in diastole, John Wiley Sons. Ltd., 30 (2014), 1199–1222. https://doi.org/10.1002/cnm.2652 doi: 10.1002/cnm.2652
    [12] L. Cai, H. Gao, X. Y. Luo, Y. H. Wang, Y. Q. Li, A mathematical model for active contraction in healthy and failing myocytes and left ventricles, Plos One, 12 (2017), e0174834. https://doi.org/10.1371/journal.pone.0174834 doi: 10.1371/journal.pone.0174834
    [13] L. C. Yann, B. Yoshua, H. Geoffrey, Deep learning, Nature, 521 (2015), 436–44. https://doi.org/10.1038/nature14539 doi: 10.1038/nature14539
    [14] G. Litjens, T. Kooi, B. E. Bejnordi, A. Setio, C. I. Sánchez, A survey on deep learning in medical image analysis, Med. Image Anal., 42 (2017), 60–88. https://doi.org/10.1016/j.media.2017.07.005 doi: 10.1016/j.media.2017.07.005
    [15] D. Shen. G. Wu, H. I. Suk, A survey on deep learning in medical image analysis, Annu. Rev. Biomed. Eng., 19 (2017), 221–248. https://doi.org/10.1146/annurev-bioeng-071516-044442 doi: 10.1146/annurev-bioeng-071516-044442
    [16] C. Case, J. Casper, G. Diamos, Deep Speech: Scaling up end-to-end speech recognition, Comput. Sci., 2014.
    [17] J. Bonnemain, L. Pegolotti, L. Liaudet, S. Deparis, Implementation and calibration of a deep neural network to predict parameters of left ventricular systolic function based on pulmonary and systemic arterial pressure signals, Front. Physiol., 11 (2022), 1086. https://doi.org/10.3389/fphys.2020.01086 doi: 10.3389/fphys.2020.01086
    [18] J. Bonnemain, M. Zeller, L. Pegolotti, S. Deparis, L. Liaudet, Deep neural network to accurately predict left ventricular systolic function under mechanical assistance, Front. Cardiovasc. Med., 8 (2021), 752088. https://doi.org/10.3389/fcvm.2021.752088 doi: 10.3389/fcvm.2021.752088
    [19] L. Cai, L. Ren, Y. H. Wang, W. X. Xie, H. Gao, Surrogate models based on machine learning methods for parameter estimation of left ventricular myocardium, R. Soc. Open Sci., 8 (2021), 201121. https://doi.org/10.1098/rsos.201121 doi: 10.1098/rsos.201121
    [20] M. L. Liu, L. Liang, W. Sun, Estimation of in vivo constitutive parameters of the aortic wall: a machine learning approach, Comput. Methods Appl. Mech. Eng., 347 (2018), 201–217. https://doi.org/10.1016/j.cma.2018.12.030 doi: 10.1016/j.cma.2018.12.030
    [21] H. Gao, W. G. Li, L. Cai, C. Berry, X. Y. Luo, Parameter estimation in a Holzapfel-Ogden law for healthy myocardium, J. Eng. Math., 95 (2015), 231–248. https://doi.org/10.1007/s10665-014-9740-3 doi: 10.1007/s10665-014-9740-3
    [22] H. Gao, A. Aderhold, K. Mangion, Changes and classification in myocardial contractile function in the left ventricle following acute myocardial infarction, J. R. Soc. Interf., 14 (2017), 20170203. https://doi.org/10.1098/rsif.2017.0203 doi: 10.1098/rsif.2017.0203
    [23] S. F. Pravdin, H. Dierckx, L. B. Katsnelson, O. Solovyova, V. S. Markhasin, A. V. Panfilov, Electrical wave propagation in an anisotropic model of the left ventricle based on analytical description of cardiac architecture, Plos One, 9 (2014), e93617. https://doi.org/10.1371/journal.pone.0093617 doi: 10.1371/journal.pone.0093617
    [24] G. A. Holzapfel, R. W. Ogden, Constitutive modelling of passive myocardium: a structurally based framework for material characterization, Philos. Trans., 367 (2009), 3445–3475. https://doi.org/10.1098/rsta.2009.0091 doi: 10.1098/rsta.2009.0091
    [25] H. M. Wang, H. Gao, X. Y. Luo, C. Berry, B. E. Griffith, R. W. Ogden, et al., Structure-based finite strain modelling of the human left ventricle in diastole, Int. J. Numer. Methods Biomed. Eng., 29 (2013), 83–103. https://doi.org/10.37059/tjosal.2013.29.1.83 doi: 10.37059/tjosal.2013.29.1.83
    [26] L. Fan, J. Yao, C. Yang, D. Xu, D. Tang, Patient-specific echo-based left ventricle models for active contraction and relaxation using different zero-load diastole and systole geometries, Int. J. Comput. Methods, 16 (2017), 1842014. https://doi.org/10.1142/S0219876218420148 doi: 10.1142/S0219876218420148
    [27] G. Viatcheslav, A high-resolution computational model of the deforming human heart, Biomech. Model. Mechanobiol., 14 (2015), 829–849. https://doi.org/10.1007/s10237-014-0639-8 doi: 10.1007/s10237-014-0639-8
    [28] Y. Xu, Simulation on the Mechanic Behaviour of Heart Based on Finite Element Methods, Ph.D thesis, Harbin Institute of Technology in Harbin, 2015.
    [29] M. D. Mckay, R. J. Conover, J. Beckmanw, A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics, 21 (1979), 239–245. https://doi.org/10.2307/1268522 doi: 10.2307/1268522
    [30] K. Stefan, Single-beat estimation of end-diastolic pressure-volume relationship: a novel method with potential for noninvasive application, Am. J. Physiol., 291 (2006), H403–H412. https://doi.org/10.1152/ajpheart.01240.2005 doi: 10.1152/ajpheart.01240.2005
    [31] G. Xavier, B. Antoine, B. Yoshua, Deep sparse rectifier neural networks, J. Mach. Learn. Res., 15 (2011), 315–323.
    [32] G. A. Montague, M. T. Tham, M. J. Willis, A. J. Morris, Predictive control of distillation columns using dynamic neural networks, IFAC Proceed. Vol., 25 (1992), 243–248. https://doi.org/10.1016/S1474-6670(17)50999-4 doi: 10.1016/S1474-6670(17)50999-4
    [33] M. D. Mckay, R. J. Conover, J. Beckmanw, Adam: A method for stochastic optimization, CoRR, 2014. https://doi.org/10.48550/arXiv.1412.6980 doi: 10.48550/arXiv.1412.6980
  • This article has been cited by:

    1. Shasha Gao, Maia Martcheva, Hongyu Miao, Libin Rong, A two-sex model of human papillomavirus infection: Vaccination strategies and a case study, 2022, 536, 00225193, 111006, 10.1016/j.jtbi.2022.111006
    2. Chaoqian Wang, Ziwei Wang, Qiuhui Pan, Ning Cai, Injurious information propagation and its global stability considering activity and normalized recovering rate, 2021, 16, 1932-6203, e0258859, 10.1371/journal.pone.0258859
    3. Qinyue Zheng, Xinwei Wang, Qiuwei Pan, Lei Wang, Optimal strategy for a dose-escalation vaccination against COVID-19 in refugee camps, 2022, 7, 2473-6988, 9288, 10.3934/math.2022515
    4. Ye Xuan Li, Hua Liu, Yu Mei Wei, Ming Ma, Gang Ma, Jing Yan Ma, Ljubisa Kocinac, Population Dynamic Study of Prey-Predator Interactions with Weak Allee Effect, Fear Effect, and Delay, 2022, 2022, 2314-4785, 1, 10.1155/2022/8095080
    5. Yen-Chang Chang, Ching-Ti Liu, A Stochastic Multi-Strain SIR Model with Two-Dose Vaccination Rate, 2022, 10, 2227-7390, 1804, 10.3390/math10111804
    6. T.A. Midhun, K. Murugesan, A study of stochastically perturbed epidemic model of HPV infection and cervical cancer in Indian female population, 2025, 228, 03784754, 431, 10.1016/j.matcom.2024.09.008
    7. Ramziya Rifhat, Zhidong Teng, Lei Wang, Ting Zeng, Liping Zhang, Kai Wang, Dynamical behavior and density function of a stochastic model of HPV infection and cervical cancer with a case study for Xinjiang, China, 2023, 360, 00160032, 7770, 10.1016/j.jfranklin.2023.06.008
    8. Hua Liu, Xiaofen Lin, Xinjie Zhu, Qibin Zhang, Yumei Wei, Gang Ma, Modeling and analysis of a human papilloma virus transmission model with impact of media, 2024, 375, 00255564, 109247, 10.1016/j.mbs.2024.109247
    9. 丽娜 王, Dynamic Analysis of a Kind of HPV Transmission Model Incorporating Media Impact and Early Screening, 2024, 13, 2324-7991, 3845, 10.12677/aam.2024.138366
    10. H.J. Alsakaji, Y.A. El-Khatib, F.A. Rihan, A. Hashish, A stochastic epidemic model with time delays and unreported cases based on Markovian switching, 2024, 6, 25889338, 234, 10.1016/j.jobb.2024.08.002
    11. Bushra Bajjah, Mahmut Modanli, Francisco R. Villatoro, Finite Difference Method for Infection Model of HPV with Cervical Cancer under Caputo Operator, 2024, 2024, 1607-887X, 1, 10.1155/2024/2580745
    12. Jingwen Xu, Guzainuer Abudurusuli, Jia Rui, Zhuoyang Li, Zeyu Zhao, Yilan Xia, Xiaohao Guo, Buasiyamu Abudunaibi, Benhua Zhao, Qiwei Guo, Jing-An Cui, Yulin Zhou, Tianmu Chen, Epidemiological characteristics and transmissibility of HPV infection: A long-term retrospective study in Hokkien Golden Triangle, China, 2013–2021, 2023, 44, 17554365, 100707, 10.1016/j.epidem.2023.100707
    13. Ramziya Rifhat, Zhidong Teng, Lei Wang, Ting Zeng, Liping Zhang, Kai Wang, Mathematical modeling analysis and simulation of human papillomavirus infection and cervical cancer in Xinjiang, China, 2023, 46, 0170-4214, 18651, 10.1002/mma.9584
    14. Ali Akgül, Nauman Ahmed, Sadiya Ali Rano, Qasem Al-Mdallal, A hybrid fractional model for cervical cancer due to human papillomavirus infection, 2025, 26662027, 101098, 10.1016/j.ijft.2025.101098
    15. Sylas Oswald, Eunice Mureithi, Berge Tsanou, Michael Chapwanya, Kijakazi Mashoto, Crispin Kahesa, MCMC-Driven mathematical modeling of the impact of HPV vaccine uptake in reducing cervical cancer, 2025, 24682276, e02633, 10.1016/j.sciaf.2025.e02633
    16. Hua Liu, Xinjie Zhu, Xiaofen Lin, Qibin Zhang, Yumei Wei, Stability analysis and optimal control of a SVICR HPV model with vaccination and cancerous delay, 2025, 1598-5865, 10.1007/s12190-024-02335-6
    17. Ramziya Rifhat, Shayidan Abuduwaili, Zhidong Teng, Kai Wang, Analysis of vaccination strategies in a heterosexual HPV transmission model with a case study in Xinjiang of China, 2025, 24680427, 10.1016/j.idm.2025.06.004
    18. Md Mehedi Hasan, Md Hamidul Islam, Junyuan Yang, Herd immunity and prevention in HPV transmission with exogenous reinfection, 2025, 20, 1932-6203, e0327233, 10.1371/journal.pone.0327233
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2580) PDF downloads(112) Cited by(0)

Figures and Tables

Figures(13)  /  Tables(9)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog