Fractional Order In Vitro Fertilization Model Real Data Analysis with Novel Application of Inequalities via Stability and Computational Techniques

1  Introduction

Infertility is an issue that, unfortunately, many couples struggle with, up to 12 percent worldwide [1]. It is defined as a disease that causes the failure of conceiving after one year of unprotected sexual intercourse, according to the International Committee for Monitoring Assisted Reproductive Technology and the World Health Organization [2]. In Vitro Fertilization (IVF) has been rapidly becoming a major Assisted Reproductive Technology (ART) to help couples facing such problems to achieve pregnancy [3–5]. The anniversary of the first “test-tube baby” birth was in 1978 on the 25th of July [6]. Different forms of infertility were addressed. These included female tubal disorders and/or inflammatory diseases, male spermatogenesis defects or insufficient sperm transport, couples experiencing long-term unexplained infertility, the presence of spermatozoa immunity in one or both spouses, and couples experiencing combined factors [7]. There are two different types of IVF cycles: regular IVF, where both the sperm and the oocyte are retrieved from the parents, and donor cycle IVF, where either the sperm or the oocyte (or both) was retrieved from a healthy, fertile donor [8]. Even so, concerns have been rising regarding the impact of factors like decreasing embryo implantation rates and increasing maternal age on IVF outcomes [9–11]. To have the first live birth following IVF and embryo transfer, there were certain challenges that had to be conquered, like ovarian stimulation, cycle monitoring, oocyte culture, sperm preparation and capacitation, etc. [12,13]. An important indicator of a successful embryo implantation after IVF is a positive beta-human chorionic gonadotropin (βHCG) [14,15].

Reproductive experts have struggled for years to predict the success of pregnancy in IVF cycles [16]. Using mathematical modeling to represent the IVF dynamics has become a very strong approach in medical research [17]. The two primary categories of mathematical models are empirical models and mechanistic models. Empirical models rely on data, which are frequently employed in medical research and are designed for prediction. Mechanistic models take into account the components of a system and how they interact. They excel in formulating, evaluating, and improving hypotheses [18].

An effective technique for analyzing an extensive diversity of real-life problems in many fields, in engineering and science, is applying the concept of a fractional calculus (FC), like Riemann-Liouville, Caputo, and Granöwald-Letnikov; all are among the theories that describe the derivative and integrated functions [19]. Since fractional calculus can model the majority of real-world problems, it has attracted a lot of attention in contemporary science and technology, in contrast to classical calculus. Numerous qualitative and numerical investigations of differential and integral equations of non-integer order have resulted from this [20–22].

Fractal-fractional operators are a new development in fractional calculus, which merge the characteristics of fractional power laws and fractal geometry to simulate even more complicated hereditary behavioral patterns [23]. This study uses the Caputo fractional derivative to preserve compatibility with the classical initial circumstances, while the fractal fractional operators show a good frontier for more improving reproductive models.

Due to its capacity to capture memory and hereditary effects that are not captured in traditional integer-order systems, fractional-order differential operators have garnered considerable interest in epidemiological modeling [24]. In biological and biomedical systems, fractional calculus has been effectively used to describe memory and inherited traits [25]. In spite of this, it has not been extensively used for predicting IVF outcomes. However, hormonal reactions, female ovarian reserve reduction, and embryo formation tracks are examples of biological procedures in reproductive health care that frequently display memory effects [26]. These impacts are not captured by typical integer order models.

This paper explores factors influencing positive βHCG rates, including the origin of oocytes (own vs. donor), embryo quality, and embryo status (fresh vs. frozen), and also evaluates the influence of the number of embryos transferred on the outcomes of the IVF cycles by analyzing the live birth rates associated with single vs. multiple embryo transfers. The main aim of this work is specifically to investigate these concerns at the Near East IVF Clinic in Northern Cyprus and, through the examination of these factors, to provide significant knowledge and understanding regarding the IVF procedures carried out at the Near East IVF Clinic. The results will have a substantial impact on providing guidance to patients and shaping the creation of efficient public policies regarding ART. Using data from IVF clinics at Near East University Hospital that describe the factors that affect the success of pregnancy, such as the male and female factors. In contrast to simply statistical methods, this model takes into account the fundamental mechanics of IVF progress, and it uses a fractional order to detect memory effects. To validate the model, actual clinical data from Near East IVF Clinic are used to help in predicting the outcomes of IVF treatment, which is dependent on numerous factors. The findings of this study may provide significant and novel insights into the basic workings of the IVF process, enabling medical professionals to tailor therapies to each patient’s specific needs and improve the chances of a successful pregnancy.

2  Materials and Methods

This paper employs a retrospective analysis of data from 253 IVF treatment cycles at Near East University Hospital between 2015 and 2023. Patient information, including male and female factors, donor details, and IVF outcomes, was obtained from clinical records approved by the Faculty of Medicine, Department of Medical Genetics, Near East University, Center of Excellence, and Mathematics Research Center. The study was conducted in accordance with relevant guidelines and regulations and received ethical approval under the number NEU/2024/121-1820. The data were anonymized, and participant identities could not be identified; therefore, informed consent was not required. The objective of the statistical analysis, conducted using IBM SPSS Statistics version 27 software, was to compare the rates of success among different groups of patients. A Chi-square test was conducted to examine the association between embryo quality categories (low, medium, high) and IVF success (positive βHCG outcome). This test determines whether embryo quality is significantly related to pregnancy outcomes. An independent samples t-test was used to compare the mean age between patients with successful IVF outcomes and those with unsuccessful outcomes, in order to assess the effect of maternal age on IVF success. A p-value is less than 0.05 was considered statistically significant.

Fig. 1 illustrates the stages that the regular IVF patients go through [27], starting with clinical and laboratory testing, setting up the IVF cycle, ordering drugs, signing paperwork, and obtaining nursing instruction, secondly, ovarian hyper stimulation from 8–12 days of drug injections to encourage the development of follicles, the source of future eggs, then, trigger injection when 3 or more dominant follicles at or over 17 mm are seen via ultrasound, to start oocyte pick-up (OPU) after 35–36 h by a quick, straightforward surgical operation carried out under anaesthesia, lasting 15 min. Now, if an oocyte is yielded in OPU, sperm is obtained, and fertilization of the oocytes via intracytoplasmic sperm injection (ICSI) technique is used, and the endometrium is supported with pure progesterone. The embryo is followed up for 2–5 days to continue its development. The last stage will be the embryo transfer, and after 12 days from the embryo transfer βHCG pregnancy test is done. The descriptive and the characteristics of the sample of the study are shown in Table 1. Different types of modalities were used in this data. Most embryo grading methods, such as the Istanbul Consensus, assign a grade of “good,” (G1), “medium,” (G2), or “poor” (G3) to embryos without employing a numerical index. The non-numerical factor-based criteria vary based on the architecture of the embryo evaluation system for each developmental stage. Therefore, using non-numerical indicators based on distinct evaluation standards makes it impossible to compare the features of groups of embryos at different developmental stages. Because of this, this type of comparison needs a numerical scoring system based on equivalent evaluation criteria [28]. From 2015 to 2023, about 253 patients of IVF treatment were included in this clinical data analysis. The data was obtained from the Near East University Hospital IVF clinic. This data, which consisted of patient files and clinical outcome records, is shown in Table 1. These are a select few variables and parameters that help highlight the factors that affect the result of IVF in Northern Cyprus.

Figure 1: Schematic illustration of the general IVF stages.

3  Formalization of Fractional Order IVF Model

It can be difficult to model IVF success because it involves many factors that affect the process itself. However, we can demonstrate this model with a particular example of a compartmental model using Caputo Fractional Calculus with six-compartment model of IVF, despite the fact that it won’t cover every situation. The six compartments of the model are:

•   I(t): the number of infertile couples that diagnosed as Polycystic Ovary Syndrome (PCOS), male factors or female factors;

•   V(t): the number of patients under going IVF treatment including ovarian stimulation, egg retrieval, sperm collection and fertilization;

•   G1(t): the number of high quality embryos;

•   G2(t): the number of medium quality embryos;

•   G3(t): the number of low quality embryos; and

•   F(t): the number of positive pregnancy results.

Here we consider some useful results of fractional calculus for model formulation and analysis for the given assumption.

Definition 1 ([25]): Assume that J is a real value function on L1([0,T],R) so the Riemann Liouville fractional order τ∈(0,1) written as:

⁢0RLDtτJ(t)=1Γ(τ)∫0t(t−μ)(τ−1)J(μ)dμ.(1)

The interval (0,∞) defines the integral on the right side.

Definition 2 ([25]): Let J be a real value function on [0,T]. then the Caputo derivative of J is:

⁢0CDtτJ(t)=1Γ(n−τ)[∫0t(t−μ)n−μ−1dndμnJ(t)(μ)dμ],(2)

where the Gamma function is represented by the symbol Γ(.)

Let n=1, then the Eq. (2) can be expressed as:

⁢0CDtτJ(t)=1Γ(n−τ)[∫0t(t−μ)−τJ′(t)(μ)dμ].(3)

Lemma 1: Let τ∗ be the Caputo systems equilibrium point, then

Dtτ0CJ(t)=J(t,τ(t)),τ∈(0,1),(4)

if and only if J(t,τ∗)=0.

Every interaction in the suggested compartmental model is directly connected to recognized clinical procedures in IVF treatment in order to improve biological interpretability. Transition from infertile couples to the second compartment represents the clinical decision and eligibility for IVF treatment. The associated rate reflects factors such as medical recommendation, patient choice, and access to fertility services, while the transitions from the compartment of IVF treatment to each of the embryo quality compartments are the result of fertilization, oocyte retrieval, ovarian stimulation, and embryo culture. Standard embryological grading systems based on morphological and developmental parameters are reflected in the distribution of embryos into low, medium, and high quality categories. The last transitions from embryo quality compartments to pregnancy outcome represent embryo transfer and implantation. The corresponding rates are interpreted as implantation probabilities, which are known to depend strongly on embryo quality, with high-quality embryos having significantly higher success rates.

A greater number of real systems are included in fractional-order systems, which are variations of conventional integer-order systems with a derivative between 0 and 1. We now extend and modify the model for examining IVF effectiveness using a set of fractional order equations inside the framework of the Caputo derivative, which is expressed as follows:

DtτCI(t)=πN−μ1I−βI,DtτCV(t)=βI−(μ2+α1+α2+α3)V,DtτCG1(t)=α1V−(γ1+μ)G1,DtτCG2(t)=α2V−(γ2+μ)G2,DtτCG3(t)=α3V−(γ3+μ)G3,DtτCF(t)=(γ1G1+γ2G2+γ3G3)11+σE−μF.(5)

The fractional-order IVF model is built using the following assumptions to guarantee biological explanation and mathematical simplicity:

1. The initial conditions are:

I(0)=I0≥0,V(0)=0≥0,G1(0)=0≥0,G2(0)=0≥0,G3(0)=0≥0,F(0)=0≥0.(6)

2. People move through the compartments only in one direction.

3. Embryos are classified into groups of poor, medium, and high quality. There is a unique and consistent implantation probability for every class.

4. The chance of obtaining a positive βHCG result is influenced separately by every embryo’s quality group.

5. A decreasing function that is added into the transition rates to the pregnancy compartment is used to describe how the chance of an effective implantation declines with maternal age.

6. A positive βHCG test, is the definition of IVF success which indicates biochemical pregnancy and does not model future phases like live birth.

7. All of the parameters are assumed to be constant and positive.

The dynamical system is depicted in Fig. 2, while Table 2 provides a description of each parameter.

Figure 2: Flow chart of dynamical system.

Where, γ1 represents the probability that an embryo of high quality produces a positive pregnancy, γ2 is the probability that an embryo of medium quality produces a positive pregnancy, and γ3 is the probability that an embryo of low quality produces a positive pregnancy.

All metrics, including embryo grading distributions, implantation rates, and treatment results, are specified in biologically meaningful ways and correlate to quantifiable clinical quantities that have been taken from an IVF clinic. The model equations’ dimensional consistency is maintained by the appropriate selection of units. Noted that the Parameters are estimated using empirical proportions from clinical data and refined through curve fitting to match observed IVF outcomes.

4  Analysis of the Model

4.1 Equilibrium Point Analysis

The equilibrium points will be checked in this section. As the system (5) is a flow system but not infectious one, then equilibrium points are:

E∗=(I∗,V∗,G1∗,G2∗,G3∗,F∗),

where

I∗=πNμ1+β,V∗=βπN(μ1+β)(μ2+α1+α2+α3),G1∗=α1βπN(γ1+μ)(μ1+β)(μ2+α1+α2+α3),G2∗=α2βπN(γ2+μ)(μ1+β)(μ2+α1+α2+α3),G3∗=α3βπN(γ3+μ)(μ1+β)(μ2+α1+α2+α3).(7)

The suggested system is considered absolutely stable. At this equilibrium point, stability means a constant flow of infertile couples, constant treatment rate with a stable embryo distribution and success rate [29], which leads to a clinical steady state.

4.2 The Reproduction Number (IVF Success Threshold)

In disease epidemiology, the basic reproduction number is a threshold quantity that establishes whether or not a disease can be controlled. When calculating this number for biological models, researchers frequently employ the next-generation matrix technique [30]. As this model is not a model of infectious disease, a definition of IVF success threshold is defined using [31,32], to be:

Rs=total successful implantation rate/loss rate

then,

Rs=α1γ1+α2γ2+α3γ3μ(1+σE).(8)

The R0 which is the basic reproduction number, indicates how likely it is for a disease to spread throughout a community. It is mathematically associated with the peak and ultimate extent of an epidemic and serves as a threshold for the stability of a disease-free equilibrium. In the suggested model, the IVF success threshold Rs will be used instead, which can be defined as the typical quantity of viable embryos produced, under optimal circumstances.

4.3 Sensitivity Analysis

Definition 3 ([31]): Sensitivity analysis is a mathematical technique used to quantify how variations in model parameters affect the behavior of the system.

For a fractional-order dynamical system, the sensitivity of a state variable y(t) with respect to a parameter r is defined as:

Sr(t)=∂y(t)∂r.

In normalized form, the sensitivity function is given by:

Srrel(t)=ry(t)⋅∂y(t)∂r.

This measure describes the relative change in the state variable due to a relative change in the parameter. In the context of fractional-order models, sensitivity analysis is particularly important due to the presence of memory effects, where parameter variations may influence the system dynamics over time.

The impact of each parameter in Rs is demonstrated by the analytical examples that follow:

∑μRs=−1<0,∑σRs=−σE(1+σE<0,∑α1Rs=α1γ1(α1γ1+α2γ2+α3γ3)>0∑α2Rs=α2γ2(α1γ1+α2γ2+α3γ3)>0∑α3Rs=α3γ3(α1γ1+α2γ2+α3γ3)>0∑γ1Rs=α1γ1(α1γ1+α2γ2+α3γ3)>0∑γ2Rs=α2γ2(α1γ1+α2γ2+α3γ3)>0∑γ3Rs=α3γ3(α1γ1+α2γ2+α3γ3)>0.(9)

A positive sensitivity index indicates that as the related parameter value increases, Rs will also increase. Conversely, a negative sensitivity index indicates that Rs will decrease as the parameter value increases. The main factors that affect the probability of IVF success can be found using any of these techniques. The study suggests that increasing parameters α1,α2,α3,γ1,γ2,andγ3 increase Rs, while increasing parameters μ, and σ, decrease it. Fig. 3 depicts simulation results for various parameter effects on Rs, which are within a feasible domain, and it clarifies the effect of each parameter on the IVF success threshold.

Figure 3: Analysis of Rs with different parameters.

4.4 Boundedness and Positivity of the System

This study looks at the conditions under which the proposed model, whose output is guaranteed to be bounded and positive, can be applied in real-world scenarios. The following are true ∀t≥0 as the classical derivatives:

I(t)≥I0e−(μ1+β)t,V(t)≥V0e−(μ2+α1+α2+α3)t,G1(t)≥G10e−(γ1+μ)t,G2(t)≥G20e−(γ2+μ)t,G3(t)≥G30e−(γ3+μ)t,F(t)≥F0e−μt.(10)

As long as the initial conditions of the Caputo fractional operator are met ∀t>0, the system (5) solution will remain positive [29], resulting in:

I(t)≥I0 eτ{−(μ1+β)tτ},V(t)≥V0 eτ{−(α3+α2+α1+μ2)tτ},G1(t)≥G10 eτ{−(γ1+μ)tτ},G2(t)≥G20 eτ{−(γ2+μ)tτ},G3(t)≥G30 eτ{−(γ3+μ)tτ},F(t)≥F0 eτ{−(μ)tτ}.(11)

4.5 The Positively Invariant Region

Theorem 1: The initial conditions that provided by the system’s solutions of (5) in the positive region R+5 are not negative.

Proof: The positive outcomes of the model (5) are:

DtτCI(t)|I=0=πN≥0,DtτCV(t)|V=0=βV≥0,DtτCG1(t)|G1=0=α1V≥0,DtτCG2(t)|G2=0=α2V≥0,DtτCG3(t)|G3=0=α3V≥0,DtτCF(t)|F=0=(γ1G1+γ2G2+γ3G3)11+σE≥0.(12)

In consequence, we have

∏~=(I,V,G1,G2,G3,F)∈R+6.(13)

The region ∏~ is positively invariant. □

4.6 The Uniqueness and the Existence of the System

This section uses the fixed-point theorem to investigate the existence and uniqueness of the proposed model. While Schauder’s fixed-point theorem ensures the model’s existence, Banach’s contraction theorem ensures the singularity of the model. If there is a solution, the appropriate approximation for it can be made using numerical approaches for 0<κ≤1. It is possible to make a system (5) more applicable by employing a fractional derivative in the Caputo sense. Let

λ1(t,I)=πN−μ1I−βI,λ2(t,V)=βI−(μ2+α1+α2+α3)V,λ3(t,G1)=α1V−(γ1+μ)G1,λ4(t,G2)=α2V−(γ2+μ)G2,λ5(t,G3)=α3V−(γ3+μ)G3,λ6(t,F)=(γ1G1+γ2G2+γ3G3)11+σE−μF.(14)

Now, apply the fractional integral and initial conditions, to get;

I(t)=I0+1Γ(τ)∫0t(t−ω)τ−1ρ1(t,I)dω,V(t)=V0+1Γ(τ)∫0t(t−ω)τ−1ρ2(t,V)dω,G1(t)=G10+1Γ(τ)∫0t(t−ω)τ−1ρ3(t,G1)dω,G2(t)=G20+1Γ(τ)∫0t(t−ω)τ−1ρ4(t,G2)dω,G3(t)=G30+1Γ(τ)∫0t(t−ω)τ−1ρ5(t,G3)dω,F(t)=F0+1Γ(τ)∫0t(t−ω)τ−1ρ6(t,F)dω.(15)

Let

Λ(t)={I(t)V(t)G1(t)G2(t)G3(t)F(t),ψ0={I0V0G10G20G30F0,ℵμ,Λ(μ)={λ1(t,I),λ2(t,V),λ3(t,G1),λ4(t,G2),λ5(t,G3),λ6(t,F).(16)

Therefore, the system (15) can be written as,

ϱ(t)=ϱ0+1Γ(τ)∫0t(t−ω)τ−1ℵ(ω),ϱ(ω)dω.(17)

The Banach space ℘[0,Q]=ε is found on a norm:

‖I,V,G1,G2,G3,F‖=maxt∈[0,Q][|I,V,G1,G2,G3,F|].(18)

Define a mapping as Y:ε→ε

Yϱ(t)=ϱ0+1Γ(τ)∫0t(t−ω)τ−1ℵ(ω),ϱ(ω)dω.(19)

Using the assumptions given bellow on a non-linear function:

(P1) A constants ηm,ηm>0 exists, so that

|ℵ(t,ϱ(t)|≤ηm|ϱ(t)|+ηm.(20)

(P2) A constant ζm>0 exists for each ϱ,ϱ¯∈ε such that

|ℵ(t,ϱ¯)−ℵ(t,ϱ¯1)|≤ζm|ϱ¯−ϱ¯1|.(21)

Theorem 2: Now, if (P1) is true, then the system (5) will contain at least one solution.

Proof: To show that Y is bounded, assume that ϕ={ϱ∈ϕ|η≥‖ϱ‖}, where η,

η≥maxt∈[0,Q]ϱ0+(ηmQτΓ(τ+1))1−(ηmQτΓ(τ+1)),(22)

is closed subset of ε.

‖Yϱ‖=maxt∈[0,Q]|ϱ0+1Γ(τ)∫0t(t−ω)τ−1ℵ(ω,ϱ(ω)dω|≤|ϱ0|+1Γ(τ)∫0t(t−ω)τ−1|ℵ(ω,ϱ(ω)|dω≤|ϱ0|+1Γ(τ)∫0t(t−ω)τ−1[ηm|ϱ(t)|+ηm]dω≤|ϱ0|+[ηm‖ϱ‖+ηm]QτΓ(τ+1)≤η.(23)

As ϱ∈ϕ⇒Y(ϕ)⊆ϕ, it proves that Y is named to be bounded. Now, let t1<t2∈[0,Q], to prove that Y is continuous.

‖Yϱ(t2)−Yϱ(t1)‖=|1Γ(τ)∫0t2(t2−ω)τ−1ℵ(ω),ϱ(ω)dω−1Γ(τ)∫0t1(t1−ω)τ−1ℵ(ω),ϱ(ω)dω|≤[ηm‖ϱ‖+ηm][t2τ−t1τ]Γ(τ+1).(24)

This proves that ‖Yϱ(t2)−Yϱ(t1)‖→0 as t2→t1, using Arzela- Ascoli theory, then the operator Y can be considered as continuous function. Thus, using Schauder’s fixed point theory on the system (5) has at least one solution. Then, the system (5) has a unique solution that obtained by applying the Banach fixed point theorem. □

Theorem 3: Let (5) gives a unique solution if the condition of (P2) are hold.

Proof: Let ϱ¯,ϱ¯1∈ε. Then

‖Y(ϱ¯)−Y(ϱ¯1)‖=maxt∈[0,Q]|1Γ(τ)∫0t(t2−ω)τ−1ℵ(ω,ϱ¯(ω)dω−1Γ(τ)∫0t(t1−ω)τ−1ℵ(ω,ϱ¯1(ω)dω|≤QτΓ(τ+1)ζm|ϱ¯−ϱ¯1|.(25)

Therefore Y can be considered as a contradiction. By using Banach fixed point theory, then the system (5) said to have a unique solution. □

Now, to show that the model actually satisfies the conditions mentioned before,

Consider the proposed IVF fractional order model written as:

DτY(t)=F(Y(t)),0<τ≤1,

where

Y(t)=(I,V,G1,G2,G3,F)T

where the function F(Y)=AY, with A being a constant matrix of biological transition rates.

Theorem 4: Assume F:Rn→Rn be a continuous function satisfying the Lipschitz condition and linear growth condition. Then, the fractional differential system shows a unique solution on a given interval [25].

Proof: To prove the conditions for the proposed model, the system can be classified to be linear with:

F(Y)=AY.

For any Y,Z∈R6, we have:

‖F(Y)−F(Z)‖=‖A(Y−Z)‖.

Using properties of matrix norms, it follows that:

‖F(Y)−F(Z)‖≤‖A‖‖Y−Z‖,

which shows that F satisfies the Lipschitz condition with Lipschitz constant L=‖A‖.

Furthermore,

‖F(Y)‖=‖AY‖≤‖A‖‖Y‖,

which Illustrates the linear growth condition.

Since both Lipschitz continuity and linear growth conditions are satisfied, The system provides a unique solution, according to classic results in fractional differential equations. □

5  Stability Analysis

Consider the fractional-order linear system:

DτY(t)=AY(t),0<τ≤1,

where

Y=(I,V,G1,G2,G3,F)T

illustrates the IVF compartments, and A is a constant matrix of biological transition rates.

Theorem 5 ([33]): The equilibrium point y∗=0 of the system is locally asymptotically stable if and only if all eigenvalues λi of the matrix A satisfy:

|arg⁡(λi)|>τπ2,∀i.

In particular, if all eigenvalues satisfy:

Re(λi)<0,

then the system is locally asymptotically stable for any 0<τ≤1.

[−μ1−β0000β−μ2−α1−α2−α30000α1−γ1−μ000α20−γ2−μ00α300−γ3−μ.](26)

After checking the eigenvalues, which are determined using Maple software 2023, it was discovered that all of them were negative real values. Based on this, it can be said that the equilibrium point is locally stable.

The Lyapunov function can be used to determine the global stability of the equilibrium points, and the crucial lemma [34] can be used to examine the global stability of the suggested system.

Lemma 2: For t≥t0, let Q∈R+ be a continuous function. Then

DtτC(K−D∗−K∗log⁡KK∗)≤(1−K∗K)DtτCK(t),K∗∈R+∀τ∈(0,1).(27)

For the Lyapunov function, {I,V,G1,G2,G3}, L<0 is the equilibrium points E∗.

Proof: Let the Volterra-type Lyapunov function:

L=L1(I−I∗−I∗log⁡II∗)+L2(V−V∗−V∗log⁡VV∗)+L3(G1−G1∗−G1∗log⁡G1G1∗)+L4(G2−G2∗−G2∗log⁡G2G2∗)+L5(G3−G3∗−G3∗log⁡G3G3∗)+L6(F−F∗−F∗log⁡FF∗),(28)

where Li,i=1,2,3,4 are the positive constants. Then Lemma 2 can be used after substituting Eq. (28) into the system.

DtτCL≤L1(I−I∗I)DtτCI+L2(V−V∗V)DtτCV+L3(G1−G1∗G1)DtτCG1+L4(G2−G2∗G2)DtτCG2+L5(G3−G3∗G3)DtτCG3+L6(F−F∗F)DtτCF.(29)

Now, substitute the derivative values in the above equation to have:

DtτCL≤L1(I−I∗I)(πN−μ1I−βI)+L2(V−V∗V)(βI−V(μ2+α1+α2+α3))+L3(G1−G1∗G1)(α1V−G1(γ1+μ))+L4(G2−G2∗G2)(α2V−G2(γ2+μ))+L5(G3−G3∗G3)(α3V−G3(γ3+μ))+L6(F−F∗F)((γ1G1+γ2G2+γ3G3)11+σE−μF).(30)

Replacing I=I−I∗,V=V−V∗,G1=G1−G1∗,G2=G2−G2∗,G3=G3−G3∗,F=F−F∗, to get;

DtτCL≤L1(I−I∗I)(πN−μ1(I−I∗)−β(I−I∗))+L2(V−V∗V)(β(I−I∗)−(V−V∗)(μ2+α1+α2+α3))+L3(G1−G1∗G1)(α1(V−V∗)−(G1−G1∗)γ1)+L4(G2−G2∗G2)(α2(V−V∗)−(G2−G2∗)γ2)+L5(G3−G3∗G3)(α3(V−V∗)−(G3−G3∗)γ3)+L6(F−F∗F)(γ1(G1−G1∗)+γ2(G2−G2∗)+γ3(G3−G3∗))11+σE−μ(F−F∗)).(31)

Let L1=L2=L3=L4=L5=1, then,

DtτCL≤(πN−πNI∗I−μ1(I−I∗)2I−β(I−I∗)2I)+(βI−βI∗−βIV∗V+βI∗V∗V−(V−V∗)2V(μ2+α1+α2+α3))+(α1V−α1V∗−α1VG1∗G1+α1V∗G1∗G1−γ1(G1−G1∗)2G1)+(α2V−α2V∗−α2VG2∗G2+α2V∗G2∗G2−γ2(G2−G2∗)2G2)+(α3V−α3V∗−α3VG3∗G3+α3V∗G3∗G3−γ3(G3−G3∗)2G3)+γ1G11+σE−γ1G1∗1+σE−γ1G1F∗(1+σE)F+γ1G1∗F∗(1+σE)F)+γ2G21+σE−γ2G2∗1+σE−γ2G2F∗(1+σE)F+γ2G2∗F∗(1+σE)F)+γ3G31+σE−γ3G3∗1+σE−γ3G3F∗(1+σE)F+γ3G3∗F∗(1+σE)F)−μ(F−F∗)2F).(32)

After simplification we get,

DtτCL≤ζ1−ζ2,(33)

where

ζ1=πN+βI+βI∗V∗V+α1V+α1V∗G1∗G1+α2V+α2V∗G2∗G2+α3V+α3V∗G3∗G3+γ1G11+σE+γ1G1∗F∗(1+σE)F+γ2G21+σE+γ2G2∗F∗(1+σE)F+γ3G31+σE+γ3G3∗F∗(1+σE)F,

and.

ζ2=πNI∗I+μ1(I−I∗)2I+β(I−I∗)2I)+βI∗+βIV∗V+(V−V∗)2V(μ2+α1+α2+α3)+α1V∗+α1VG1∗G1+γ1(G1−G1∗)2G1+α2V∗+α2VG2∗G2+γ2(G2−G2∗)2G2)+α3V∗+α3VG3∗G3+γ3(G3−G3∗)2G3+γ1G1∗1+σE+γ1G1F∗(1+σE)F+γ2G2∗1+σE+γ2G2F∗(1+σE)F+γ3G3∗1+σE+γ3G3F∗(1+σE)F+μ(F−F∗)2F.

If

ζ1<ζ2   ⇒  DtτCL<0.

But, if I=I∗, V=V∗, G1=G1∗, G2=G2∗, G3=G3∗, F=F∗ then

ζ1−ζ2=0,    ⇒   DtτCL=0.

Therefore, we can conclude that if ζ1<ζ2, the equilibrium points E∗ are globally asymptotically stable in the invariant region. □

6  Numerical Algorithm for Computational Analysis

Caputo’s derivative can be used to model power-law processes, as is well known from the literature. In order to incorporate power law effects into our model, we evaluated the dynamics seen in fractional calculus using the temporal derivative in conjunction with the Caputo derivative. A numerical method based on Newton polynomial interpolation can then be used to discretize the problem at hand [29,35]. One way to express the system (5) is:

DtτCI(t)=λ1(t,Θ(t)),DtτCV(t)=λ2(t,Θ(t)),DtτCG1(t)=λ3(t,Θ(t)),DtτCG2(t)=λ4(t,Θ(t)),DtτCG3(t)=λ5(t,Θ(t)),DtτCF(t)=λ6(t,Θ(t)).(34)

where Θ(t)=(I(t),V(t),G1(t),G2(t),G3(t),F(t)).

Let the general description of (34) be

{Dtτ0CJ(t)=F(t,J(t)),τ∈(0,1],t∈[0,Q],J(0)=J0.(35)

Caputo Integral Representation

The Caputo derivative can be written as:

J(t)=J0+1Γ(τ)∫0t(t−ζ)τ−1F(J(ζ))dζ.

Time Discretization

Let the interval [0,T] be divided into uniform steps:

tκ=κh,κ=0,1,2,…,N,

where h is the step size.

Then:

Jκ+1=J0+1Γ(τ)∑n=0κ∫tntn+1(tκ+1−ζ)τ−1F(J(ζ))dζ.

Newton Polynomial Interpolation

On each interval [tn,tn+1], we approximate F(j(ζ)) using the first-order Newton interpolation polynomial:

F(J(ζ))≈F(Jn)+ζ−tnh(F(Jn+1)−F(Jn)).

Substitution into the Integral

Substituting the interpolation into the integral yields:

∫tntn+1(tκ+1−ζ)τ−1F(J(ζ))dζ=F(Jn)Aκ,n+F(Jn+1)Bκ,n,

where the coefficients are given by:

Aκ,n=∫tntn+1(tκ+1−ζ)τ−1(1−ζ−tnh)dζ,

Bκ,n=∫tntn+1(tκ+1−ζ)τ−1ζ−tnhdζ.

After simplification, these coefficients can be written as:

Aκ,n=hτΓ(τ+2)[(κ+1−n)τ+1−(κ−n)τ+1],

Bκ,n=hτΓ(τ+2)[(κ+1−n)τ+1−(κ−n)τ+1−τ(κ−n)τ].

Final Discrete Scheme

Thus, the numerical solution is given by:

Jκ+1=J0+1Γ(τ)∑n=0κ[Aκ,nF(Jn)+Bκ,nF(Jn+1)].

Component-wise System

Applying this to the IVF compartments:

Iκ+1=I0+1Γ(τ)∑n=0κ[Aκ,n(πN−μ1In−βIn)+Bκ,n(πN−μ1In+1−βIn+1)]