Global Piecewise Analysis of HIV Model with Bi-Infectious Categories under Ordinary Derivative and Non-Singular Operator with Neural Network Approach

1  Introduction

Viruses of HIV lie and transmit both in classical and modern world contents having various rates of infection. Among the sources of spreading of the said diseases from one man to another are through sexual meetings, used needles re-usages, drug usage, inherited through the mother, and by un-tested blood donations or transfusions. For stabilizing HIV infection transmission, there are various approaches including need-based HIV testing, usage of condoms during sex, circumcisions, usage of vaginal microbicides and antire-troviral (ARV) drug. To date, no cure or vaccination for HIV has been successfully discovered [1].

HIV infection is among the badly killing diseases, causing thousands of deaths all over the globe. Nearly thirty-eight million humans all over the world suffered from HIV in the year of 2019. After that, the virus has spread to each country of the globe. According to the statistics of the Health Ministry about 21,500 men were suffering from HIV in the African country of Morocco [1]. Among them, 6000(22%) were ignorant of their infection. In 2019 about 850 novel HIV testing in Morocco, 33% were aged 15–25. This analysis also gives 300 deaths data. In more detail, 67% of novel cases occur in highly risky chain populations. About 70.7% of women are tested positively due to their husbands. Statistical data implies a low percentage of HIV infection in the country of Morocco in common lie population, (0.08%),1.7% are tested females working in different fields, 5.9% on mean are the homo-sexuals, 7.1% on mean are those taking drugs by injections. Some territories of the country numbered about 65% of all infections, namely Sous-Massa (26%), Marrakech-Safi (20%) and Casablanca-Settat (20%) [1].

In dealing with such infection, one of the significant tools is mathematical modeling whose benefits are to provide an esteemed prediction and feedback of any infection transmission whose viruses are contagious and invisible. Mathematical models of HIV are constructed by the scientific approval, given by analyzers, physicians and pharmacists along with the limit of spreading of the virus. So different HIV systems have been investigated by many authors [2–5]. Researchers have examined various models by applying different methodologies and explored the behaviors of each class and comparisons to real data [6–9]. Fractional calculus can be a useful technique in HIV disease modeling to capture the complex dynamics of the illness, particularly in areas where standard models are inadequate. It does, however, also provide new difficulties about the interpretation and complexity of mathematics [10–13]. Silva et al. [12] have considered an epidemiological system of HIV/AIDS transmissions including Pre-exposed Prophylaxis. In the same fashion, Li et al. [5] made a susceptible exposure in the Latent stage infection (SEI) level to draw the evolution of HIV. One of the Vivo deterministic problems has been discussed by Ngina et al. [7]. Shirazian et al. [11] have considered a mathematical method applied to express mathematically the procedure for medical tests and seeing the infection of HIV/AIDS.

The scholars have proposed a continuous system of five agents SI1I2QR for discussion of the interaction among the populations having HIV/AIDS. Here S,I1,I2,Q,R are the cases of susceptible peoples, Infectious peoples stage 1, infectious peoples stage 2, Infectious peoples who are hospitalized and recoverable individuals, respectively [1]. The said model is given below:

dS(t)dt=Δ−uS−(β1I1+β2I2)SdI1(t)dt=β1I1S+δ2I2−(δ1+γ1+u+σ1)I1dI2(t)dt=β2I2S+δ4Q+δ1I1−(γ2+δ2+δ3+u+σ2)I2dQ(t)dt=δ3I2−(γ3+δ4+u+σ3)QdR(t)dt=σ1I1+σ2I2+σ3Q−uRS(0)=S0,I1(0)=I10,I2(0)=I20,Q(0)=Q0,R(0)=R0≥0.(1)

The used parameters are defined in Table 1 [1] along with numerical values.

Many sorts of derivatives have been developed in multiple attempts to understand cross-over dynamics and overcome discontinuities. These comprise fractal-fractional derivatives, derivatives of fractional order whose kernels show singularity and non-singularity, and other special cases of derivatives [14–17]. These methods, which offer sophisticated instruments for capturing the intricacies present in such dynamic processes, have all been designed to address particular facets of the behavior of the system [18–22]. These various derivative formulations enable a more sophisticated analysis of cross-over dynamics, allowing for the exploration of system transitions between multiple states or phases. Derivatives make it possible to model complicated, real-world processes more accurately and adaptably by taking into account both continuity and discontinuity [23–26]. This is particularly important in the domain of disease modeling, where a precise understanding of the dynamics of transitions, such as the HIV’s transition from latency to active infection, can result in more accurate forecasts and more potent interventions [27,28]. For analysis of the randomness or probabilistic in the scheme of stochastic equations having more real findings but up to now the cross-over dynamical behavior has not been investigated. Such behaviors are found in most infection models, like the heat flow, fluid dynamics, and many of the complex geometrical problems [22,29,30]. In fractional derivatives, the exponential and Mittag-Leffler mapping are not used to compute the time of crossover dynamics. Hence to treat such models, one of the novel schemes of piece-wise derivatives and antiderivatives has been formulated in [31]. A comprehensive framework incorporating global and classical piecewise derivatives was built by the concept’s authors. These derivatives were developed to model complicated systems, namely the dynamics of HIV infection. Keeping in view the said aspects, we will also discuss the said problem for at least one solution, uniqueness of solution, numerical solution, and stability analysis in the sense of classical and Atangana-Baleanu (AB) piecewise derivative. Further, such kind of operators are different and better than the other fractional operators. It removes the discontinuity on the whole interval by converting it into a sub-interval, it describes the crossover dynamics, AB operator has the non-singular kernel which removes the singularity of the given domain, it quickly gains the stability of any dynamical system on small fractional order, is more generalized than classical and other fractional operators. Moreover, the neural network (NN) approach is utilized and find different data sets for the considered model and the comparison of fractional and NN as well. The Eq. (1) is expressed in piece-wise derivative in the sense of classical and non-singular kernel operator as under.

D0ψ0POAB(S)(t)=Δ−uS−(β1I1+β2I2)SD0ψ0POAB(I1)(t)=β1I1S+δ2I2−(δ1+γ1+u+σ1)I1D0ψ0POAB(I2)(t)=β2I2S+δ4Q+δ1I1−(γ2+δ2+δ3+u+σ2)I2D0ψ0POAB(Q)(t)=δ3I2−(γ3+δ4+u+σ3)QD0ψ0POAB(R)(t)=σ1I1+σ2I2+σ3Q−uR,t∈[0,T], 0≤t≤t1, t1≤t≤t2,(2)

where POAB is for piece-wise ordinary and AB derivative in two sub-intervals of [0,T]. In more simplest form we can write Eq. (2) under

D0ψ0POAB(S(t))={Dt0O(S(t))=ddtF1(S,t),0<t≤t1,Dtψ0AB(S(t))=ABF1(S,t)t1<t≤t2,D0ψ0POAB(I1(t))={Dt0O(I1(t))=ddtF2(I1,t),0<t≤t1,Dtψ0AB(I1(t))=ABF2(I1,t),t1<t≤t2,D0ψ0POAB(t))={Dt0O(I2(t))=ddtF3(I2,t),0<t≤t1,Dtψ0AB(I2(t))=ABF3(I2,t),t1<t≤t2,D0ψ0POAB(Q(t))={Dt0O(Q(t))=ddtF4(Q,t),0<t≤t1,Dtψ0AB(Q(t))=ABF4(Q,t),t1<t≤t2,D0ψ0POAB(R(t))={Dt0O(R(t))=ddtF5(R,t),0<t≤t1,Dtψ0AB(R(t))=ABF5(R,t),t1<t≤t2,(3)

where Dt0O and Dtψ0AB are ordinary and AB fractional derivative respectively and Fi from i=1,2,3,4,5 is for left side of Eq. (2).

1.1 Basis Properties of the Model

For any epidemiological problem, most of the authors assume that the total size of population N(t) is constant, in the rest of the paper, i.e.,

S+I1+I2+Q+R=N(t)t≥0

The transfer diagram may also be seen in [1].

Lemma 1.1. The set

Ξ={S,I1,I2,Q,R}>0∈R5andS+I1+I2+Q+R<Δu

is bounded in positive feasible region.

Proof. Adding all the five agents of Eq. (2) as follows:

S+I1+I2+Q+R=N(t)t≥0,

or

Dtψ0POAB(N(t))=Δ−u(S+I1+I2+Q+R)−γ1I1−γ2I2−γ3Q≤Δ−uN(t)

Dtψ0POAB(N(t))+uN(t)≤Δ.

On application of piece-wise integration we get

It0ON(t)≤{N(0)e−ut+Δu(1−e−ut),0<t≤t1,N(t1)+1−ψABψN(t)+ψABψΓψ∫t1t(t−τ)ψ−1N(τ)d(τ)t1<t≤t2.

In the first interval if t1 is large then N(t)≤Δu∀t≥0 and in the second interval if t2→∞, then N(t)≤Δu,∀t1≥0. Hence for both cases the system is bounded in the feasible region. ■

Lemma 1.2. If the initial values of all the agent of the (2) is non-negative, i.e., S(0)≥0,I1(0)≥0,I2(0)≥0,Q(t)≥0,R(0)≥0 then the solution of problem (2) for each quantity will be positive for t1>0 and t2>0.

Proof. The proof can be seen in [1]. ■

The free equilibrium point for the model (2) is E0=(S0,I10,I20,Q0,R0)=(Δu,0,0,0,0).

The Basic reproduction number for (2) [1] is

R0=ζ2((β1P2+β2P1)+(β1P2−β2P1)+4β1β2γ1γ2),

where P1=δ1+γ1+u+σ1,P2=δ2+δ3+γ2+u+σ2 and ζ=1P1P2−γ1γ2. Further if R0<1 on equilibrium point then problem (2) is locally asymptotically stable and unstable if R0>1.

R0<1 means that any infected population may transfer the disease to less than one new infected individual, which implies that HIV infection does not spread in the population.

R0>1 means that any infected population may transfer the disease to more than one new infected individual, which implies that HIV infection spreads in the population.

2  Basic Results

Now, in the next part of the article, we provide a few background definitions of the Caputo operator of derivative and integration along with the piecewise derivative concept.

Definition 2.1. The AB derivative of a function Ψ(t) under the condition Ψ(t)∈ℋ1(0,τ) is defined as follows:

Dtψ0AB(Ψ(t))=AB(ψ)1−ψ∫0tddτΨ(τ)Eψ[−ψ1−ψ(t−τ)ψ]dτ.(4)

While the integral can be written as

Dtψ0ABΨ(t)=1−ψAB(ψ)Ψ(t)+ψAB(ψ)Γ(ψ)∫0t(t−τ)ψ−1Ψ(τ)dτ.(5)

Definition 2.2. Choosing ω(t) be differentiable and g(t) is increasing function then classical piece-wise derivative [31] can be given as

It0OΨ(t)={Ψ(t),0<t≤t1,Ψ′(t)g′(t)t1<t≤t2=T,

and the integration is

It0OΨ(t)={∫0tΨ(τ)dτ,0<t≤t1,∫t1tΨ(τ)g′(τ)d(τ) t1<t≤t2,

here It0OΨ(t) and Dt0OΨ(t) are for ordinary derivative and integration for 0<t≤t1 and global derivative integration for t1<t≤t2

Definition 2.3. Let Ψ(t) be differentiable then ordinary and fractional piece-wise derivative [31] is defined as

Dtψ0OΨ(t)={Ψ′(t),0<t≤t1,Dtψ0ABΨ(t) t1<t≤t2,

while the integration can be

It0OΨ(t)={∫0tΨ(τ)dτ,0<t≤t1,1−ψABψΨ(t)+ψABψΓψ∫t1t(t−τ)ψ−1Ψ(τ)d(τ) t1<t≤t2,

here Dtψ0OΨ(t) and It0OΨ(t) are for ordinary derivative and integration for 0<t≤t1 and AB fractional derivative integration for t1<t≤T−2.

Lemma 2.1. The solution of piece-wise derivable equation

Dtψ0OABΨ(t)=F(t,Ψ(t)),0<r≤1

is

Ψ(t)={Ψ0+∫0tΨ(τ)dτ, 0<t≤t1Ψ(t1)+1−ψABψΨ(t)+ψABψΓψ∫t1t(t−τ)ψ−1Ψ(τ)d(τ), t1<t≤t2.

3  Existence and Uniqueness

This section will evaluate whether the piece-wise derivable problem under consideration has a unique solution and whether it exists or not. Based on the information provided in Lemma 2.1, we write the system (2) for this, and we further describe it as follows:

Dtψ0POABΨ(t)=F(t,Ψ),0<ψ≤1

is

Ψ(t)={Ψ0+∫0tF(s,Ψ(s))ds, 0<t≤t1Ψ(t1)+1−ψABψΨ(t)+ψABψΓψ∫t1t(t−s)ψ−1Ψ(s)ds, t1<t≤t2,(6)

Ψ(t)={S(t)I1(t)I2(t)Q(t)R(t), Ψ0={S0I10I20Q0R0, Ψt1={S(t1)I1(t1)I2(t1)Q(t1)R(t1), F(t,Ψ(t))={F1={ddtF1(S,t)F1AB(S,t)F2={ddtF2(I1,t)F2AB(I1,t)F3={ddtF3(I2,t)F3AB(I2,t)F4={ddtF4(Q,t)F4AB(Q,t)F5={ddtF5(R,t)F5AB(R,t)(7)

Consider ∞>t2≥t>t1>0 with closed norm space as G1=C[0,T] having norm

‖Ψ‖=maxt∈[0,T]|Ψ(t)|.

For the required result, we take the growth condition on the non-linear operator as

(A1) ∃ LΨ>0; ∀ G,Ψ¯∈E we have

|F(t,Ψ)−F(t,Ψ¯)|≤LF|Ψ−Ψ¯|,

(A2) ∃ CF>0 & MF>0,;

|F(t,Ψ(t))|≤CF|Ψ|+MF.

Theorem 3.1. Let F be a piece-wise continuous mapping on sub-interval 0<t≤t1 and t1<t≤t2 on [0,T], also satisfying (A2), then piece-wise derivative (3) has at least one solution on each sub-interval.

Proof. Let us assume a closed subset in both sub-intervals of 0,T as E of B by using the Schauder fixed point theorem

E={Ψ∈B:|Ψ‖≤R1,2, R>0},

Next, let us suppose an operator T:E→E and applying (6) as

T(Ψ)={Ψ0+∫0tF(τ,Ψ(τ))dτ, 0<t≤t1Ψ(t1)+1Γ(ψ)∫0t1F(τ,Ψ(τ))(t−τ)ψ−1dτ, t1<t≤t2Ψ(t1)+1−ψAB(ψ)F(t,Ψ(t))+ψAB(ψ)Γ(ψ)∫t1t(t−τ)ψ−1F(τΨ(τ))d(τ),  t1<t≤t2.(8)

On some Ψ∈E, we follows as

|T(Ψ)(t)|≤{|Ψ0|+∫0t1|F(τ,Ψ(τ))|dτ,|Ψ(t1)|+|1−ψAB(ψ)F(t,Ψ(t))+ψAB(ψ)Γ(ψ)∫t1t(t−τ)ψ−1|F(τΨ(τ))d(τ),≤{|Ψ0|+∫0t1[CF|Ψ|+MG]dτ,|Ψ(t1)|+1−ψAB(ψ)[CF|Ψ|+MF+ψAB(ψ)Γ(ψ)∫t1t(t−τ)ψ−1[CF|Ψ|+MF]dτ,≤{|Ψ0|+t1[CF|Ψ|+MF]=R1, 0<t≤t1,|Ψt1|+1−ψAB(ψ)[CF|Ψ|+MF]+ψ(T−T)ψAB(ψ)Γψ+1[CF|Ψ|+MF]dτ=R2, t1<t≤t2,≤{R1, 0<t≤t1,R2, t1<t≤t2.

The final equation suggests that T(E)⊂E since Ψ∈E. It follows that T is closed and complete. Moving forward, we will assume perfect continuity. Specifically, we will use tm<tn∈[0,t1] as the first interval of the classical derivative notion.

|T(Ψ)(tn)−T(Ψ)(tm)|=|∫0tnF(τ,Ψ(τ))dτ−∫0tmF(τ,Ψ(τ))dτ|≤∫0tn|F(τ,Ψ(τ))|dτ−∫0tm|F(τ,Ψ(τ))|dτ≤[∫0tn(CF|Ψ|+MF)−∫0tm(CF|Ψ|+MF)≤(CFΨ+MF)[tn−tm].(9)

Next from (9), we obtain tm→tn, then

|T(Ψ)(tn)−T(Ψ)(tm)|→0,as tm→tn.

Thus, in the interval [0,t1], T is equi-continuous. Next, we take the other AB sense interval ti,tj∈[t1,T] as

|𝒯(Ψ)(tj)−𝒯(Ψ)(ti)|=|1−ψAB(ψ)F(t,Ψ(t))+ψAB(ψ)Γ(ψ)∫t1tj(tj−τ)ψ−1F(τ,Ψ(τ))dτ,−1−ψAB(ψ)F(t,Ψ(t))+(ψ)AB(ψ)Γ(ψ)∫t1ti(ti−τ)ψ−1F(τ,Ψ(τ))dτ|≤ψAB(ψ)Γ(ψ)∫t1ti[(ti−τ)ψ−1−(tj−τ)ψ−1]|F(τ,Ψ(τ))|dτ+ψAB(ψ)Γ(ψ)∫titj(tj−τ)ψ−1|F(τ,Ψ(τ))|dτ≤ψAB(ψ)Γ(ψ)[∫t1ti[(ti−τ)ψ−1−(tj−ψ)ψ−1]dτ+∫titj(tj−τ)ψ−1dτ](CF|Ψ|+MF)≤ψ(CHΨ+MH)AB(τ)Γ(ψ+1)[tjψ−tiψ+2(tj−ti)ψ].(10)

Further from (10) we get tm→tn, then

|T(Ψ)(tn)−T(Ψ)(tm)|→0,as tm→tn.

Therefore, in the interval [t1,t2], T is equi-continuous. It follows that T is an equi-continuous mapping. The Arzelá-Ascoli Theorem shows that operator T is limited, uniform, and continuous. There is at least one solution to the piecewise derivable problem (3) on each subinterval, according to Schauder’s fixed point theorem. ■

Theorem 3.2. If T is a contraction operator, then the suggested piece-wise model has a unique root with condition (A1).

Proof. As previously, given a piece-wise continuous mapping T:E→E, let Ψ and Ψ¯∈E on [0,t1] be understood as

‖T(Ψ)−T(Ψ¯)‖=maxt∈[0,t1]|∫0t1F(τ,Ψ(τ))dτ−∫0t1F(τ,Ψ¯(τ))dτ|≤t1LF‖Ψ−Ψ¯‖.(11)

From (11), we have

‖T(Ψ)−T(Ψ¯)‖≤t1LF‖Ψ−Ψ¯‖.(12)

It follows that T is contraction. In summary, the problem under consideration has a unique solution inside the specified subinterval according to the outcome of the Banach contraction theorem. Additionally, in the context of the AB derivative, for the second sub-interval t∈[t1,t2] as

‖T(Ψ)−T(Ψ¯)‖≤1−ψAB(ψ)LF‖Ψ−Ψ¯‖+ψ(T−Tψ)AB(ψ)Γ(ψ+1)LF‖Ψ−Ψ¯‖.(13)

or

‖T(Ψ)−T(Ψ¯)‖≤LF[1−ψAB(ψ)+ψ(T−T)ψAB(ψ)Γ(ψ+1)]‖Ψ−Ψ¯‖.(14)

The presented problem, therefore, has a unique solution in the second sub-interval according to the Banach contraction theorem since T is contraction. Consequently, each sub-interval of the piecewise derivable model has a unique solution according to Eqs. (12) and (14). ■

Ulam Hyers Stability

To obtain the Ulam-Hyers stability for the proposed model, we need to prove that a small perturbed term in the starting condition or the parameters of the system implies small perturbations in the solution of the system. It can be achieved by deriving that the operator T is Lipschitz continuous w.r.t the starting conditions or any of the parameters of the system.

Definition 3.1. The proposed model (2) is called Ulam-Hyers stable, if for all N>0, then the following inequality holds true.

|DtψPOABCΨ(t)−F(t,Ψ(t))|<N,t∈T.(15)

Let we have a unique solution Ψ¯∈Y that is constant A>0,

||Ψ−Ψ¯||Y≤AN,forall,t∈T.(16)

Next, if we chose an increasing operator U:[0,∞)→R+, the inequality above can be written as under:

||Ψ−Ψ¯||Z≤AU(N),for each,t∈T.

If U(0)=0, then the final solution is general Ulam-Hyers (G-H-U) stable.

Remark 1. Assume that an operator U∈C(𝒯) is independent of Ψ∈W and holds U(0)=0, then

|U(t)|≤N,t∈TDtψPOABCΨ(t)=F(t,Ψ(t))+U(t),t∈T.

Lemma 3.1. Let the mapping

Dψt0POABCΨ(t)=F(t,Ψ(t)),0<ψ≤1.(17)

The solution of (17) is

Ψ(t)={Ψ0+∫0tF(z,Ψ(z))dz, 0<t≤t1Ψ(t1)+1−ψAB(ψ)F(t,Ψ(t))+ψAB(ψ)Γ(ψ)∫t1t(t−z)ψ−1F(z,Ψ(z))dz,t1<t≤T,(18)

||F(Ψ)−F(Ψ¯)||≤{Ψ0+t1[CF|Ψ|+MF]N,t∈(0,t1][(1−ψ)Γ(ψ)+(t2ψ)AB(ψ)Γ(ψ)]N=ΛN,t∈[t1,t).(19)

Theorem 3.3. Lemma (3.1) implies that if LFtψΓ(ψ)<1, then the solution to model (3) is Hyers-Ulam stable as well as G-H-U stable.

Proof. If Ψ∈𝒲 is a solution of (3) and Ψ¯∈𝒲 is also a unique solution of (3), then we have

Case: 1 for t∈(0,t], we have

||Ψ−Ψ¯||=supt∈(0,t]⁡|Ψ−(Ψ∘+∫0t1F(z,Ψ(z))dz)|≤supt∈(0,t)⁡|Ψ−Ψ∘|+t1[CF|Ψ|+MF]N(20)

On more calculation

||Ψ−Ψ¯||≤(t1[CF|Ψ|+MF])N(21)

Case: 2

||Ψ−Ψ¯||≤supt∈[t1,t)⁡|Ψ−[Ψ(t1)+1−ψAB(ψ)[F(t,Ψ(t)),]+ψAB(ψ)Γ(ψ)[∫t1t(t−z)ψ−1F(z,Ψ¯(z))dz]]|+supt∈[t1,t)⁡1−ψAB(ψ)|F(t,Ψ(t))−F(t,Ψ¯(t),)|+supt∈[t1,t)⁡ψAB(ψ)Γ(ψ)∫t1t(t−z)ψ−1|F(z,Ψ(z))−F(z,Ψ¯(z))|dz.

Using Λ=[(1−ψ)Γ(ψ)+t2ψAB(ψ)Γ(ψ)] and further calculation we have

||Ψ−Ψ¯||𝒲≤Λ(22)

N+ΛLF||Ψ−Ψ¯||𝒲(23)

We have

||Ψ−Ψ¯||𝒲≤(Λ1−ΛLF)β||Ψ−Ψ¯||𝒲.

We use

A=max{(t1Γ(ψ+1)1−LFt1Γ(ψ+1)),Λ1−ΛLF1−MF}

Now, by Eqs. (21) and (22), we have

||Ψ−Ψ¯||𝒲≤AN,ateveryt∈[t1,t).

Thus, we finally say that the solution of (3) is Hyers-Ulam stable. Furthermore, if we write N with U(N) in (24), we get

||Ψ−Ψ¯||𝒲≤AU(N),ateacht∈[t1,t).

So, we can say that the solution to the consider model (3) is G-H-U stable based on the fact that U(0)=0.

4  Numerical Procedure

In this part, we will now construct a numerical approach for the chosen piece-wise derivable model (2). For each of the two [0,T] sub-intervals, we will develop a numerical scheme in both ordinary and AB derivative forms. The piece-wise model’s numerical scheme will resemble the integer order scheme seen in [31]. Applying Eq. (2) for both ordinary and AB format using the piece-wise integral as follows:

S(t)={S0+∫0t1F1(τ)dτ,0<t≤t1,S(t1)+1−ψAB(ψ)F1AB(t)+ψAB(ψ)Γ(ψ)∫t1t(t−τ)τ−1F1AB(τ)dτ, t1<t≤t2,I1(t))={I10+∫0t1F2(τ)dτ,0<t≤t1,I1(t1)+1−ψAB(ψ)F2AB(t)+ψAB(ψ)Γ(ψ)∫t1t(t−τ)τ−1F2AB(τ)dτ, t1<t≤t2,I2(t))={I20+∫0t1F3(τ)dτ,0<t≤t1,I2(t1)+1−ψAB(ψ)F3AB(t)+ψAB(ψ)Γ(ψ)∫t1t(t−τ)τ−1F3AB(τ)dτ, t1<t≤t2,Q(t))={Q0+∫0t1F4(τ)dτ,0<t≤t1,Q(t1)+1−ψAB(ψ)F4AB(t)+ψAB(ψ)Γ(ψ)∫t1t(t−τ)τ−1F4AB(τ)dτ, t1<t≤t2,R(t))={R0+∫0t1F5(τ)dτ,0<t≤t1,R(t1)+1−ψAB(ψ)F5AB(t)+ψAB(ψ)Γ(ψ)∫t1t(t−τ)τ−1F5AB(τ)dτ, t1<t≤t2.(24)

Subsequently, we will demonstrate the methodology for model (24)’s first equation, and the remaining agents will follow the same process.

On t=tn+1

S(tn+1))={S0+∫0t1F1(S,I1,I2,Q,R,τ)dτ,0<t≤t1,S(t1)+1−ψAB(ψ)F1AB(S,I1,I2,Q,R,tn)+ψAB(ψ)Γ(ψ)∫t1tn+1(t−τ)ψ−1F1AB(τ)dτ,t1<t≤t2.(25)