
Citation: P. Mora-Raimundo, M. Manzano, M. Vallet-Regí. Nanoparticles for the treatment of osteoporosis[J]. AIMS Bioengineering, 2017, 4(2): 259-274. doi: 10.3934/bioeng.2017.2.259
[1] | Sumati Kumari Panda, Abdon Atangana, Juan J. Nieto . Correction: New insights on novel coronavirus 2019-nCoV/SARS-CoV-2 modelling in the aspect of fractional derivatives and fixed points. Mathematical Biosciences and Engineering, 2022, 19(2): 1588-1590. doi: 10.3934/mbe.2022073 |
[2] | Tahir Khan, Roman Ullah, Gul Zaman, Jehad Alzabut . A mathematical model for the dynamics of SARS-CoV-2 virus using the Caputo-Fabrizio operator. Mathematical Biosciences and Engineering, 2021, 18(5): 6095-6116. doi: 10.3934/mbe.2021305 |
[3] | A. D. Al Agha, A. M. Elaiw . Global dynamics of SARS-CoV-2/malaria model with antibody immune response. Mathematical Biosciences and Engineering, 2022, 19(8): 8380-8410. doi: 10.3934/mbe.2022390 |
[4] | Somayeh Fouladi, Mohammad Kohandel, Brydon Eastman . A comparison and calibration of integer and fractional-order models of COVID-19 with stratified public response. Mathematical Biosciences and Engineering, 2022, 19(12): 12792-12813. doi: 10.3934/mbe.2022597 |
[5] | Rahat Zarin, Usa Wannasingha Humphries, Amir Khan, Aeshah A. Raezah . Computational modeling of fractional COVID-19 model by Haar wavelet collocation Methods with real data. Mathematical Biosciences and Engineering, 2023, 20(6): 11281-11312. doi: 10.3934/mbe.2023500 |
[6] | A. M. Elaiw, Raghad S. Alsulami, A. D. Hobiny . Global dynamics of IAV/SARS-CoV-2 coinfection model with eclipse phase and antibody immunity. Mathematical Biosciences and Engineering, 2023, 20(2): 3873-3917. doi: 10.3934/mbe.2023182 |
[7] | Biplab Dhar, Praveen Kumar Gupta, Mohammad Sajid . Solution of a dynamical memory effect COVID-19 infection system with leaky vaccination efficacy by non-singular kernel fractional derivatives. Mathematical Biosciences and Engineering, 2022, 19(5): 4341-4367. doi: 10.3934/mbe.2022201 |
[8] | Adnan Sami, Amir Ali, Ramsha Shafqat, Nuttapol Pakkaranang, Mati ur Rahmamn . Analysis of food chain mathematical model under fractal fractional Caputo derivative. Mathematical Biosciences and Engineering, 2023, 20(2): 2094-2109. doi: 10.3934/mbe.2023097 |
[9] | M. Botros, E. A. A. Ziada, I. L. EL-Kalla . Semi-analytic solutions of nonlinear multidimensional fractional differential equations. Mathematical Biosciences and Engineering, 2022, 19(12): 13306-13320. doi: 10.3934/mbe.2022623 |
[10] | Noura Laksaci, Ahmed Boudaoui, Seham Mahyoub Al-Mekhlafi, Abdon Atangana . Mathematical analysis and numerical simulation for fractal-fractional cancer model. Mathematical Biosciences and Engineering, 2023, 20(10): 18083-18103. doi: 10.3934/mbe.2023803 |
Differential equations, especially those possessing non-negative solutions, play an important role in many areas of sciences such as biology [1,2,3,4,5], chemistry [6,7,8], epidemiology [9,10,11] or population dynamics [12,13,14] to name a few of these fields. Since we focus on examining a time-discrete, mathematical model of epidemiology in our work, we limit ourselves to a concise review of literature in this area.
In this article, we want to consider the time-continuous non-autonomous susceptible-infected-recovered (SIR) model
{S′(t)=−α(t)⋅S(t)⋅I(t)N,I′(t)=α(t)⋅S(t)⋅I(t)N−β(t)⋅I(t),R′(t)=β(t)⋅I(t),S(0)=S0,I(0)=I0,R(0)=R0 | (1.1) |
in epidemiology thoroughly investigated by Wacker and Schlüter in 2020 [15]. We provide further assumptions regarding this time-continuous model in Section 2. In 1927, one of the most famous papers in mathematical epidemiology, written by Kermack and McKendrick, was published [16]. Previously, Ross and Hudson proposed important foundations for epidemiological models [17,18]. Additional well-written books and review articles on mathematical epidemiology are authored by Hethcote [19] and Brauer and Castillo-Chávez [20]. Since the establishment of the classical SIR model, many modifications have been suggested and, due to the on-going COVID-19 pandemic, many articles have been published with respect to this topic and we can only summarize a brief selection of these works.
To this day, there are works which provide closed solution formulas for simple models in epidemiology [21,22,23]. However, more articles consider numerical approximations due to appearing nonlinearities. Especially during the COVID-19 pandemic, different approaches have been published [15,24]. In 2022, Chen and co-authors suggested a susceptible-exposed-infected-unreported-recovered (SEIUR) model for the spread of COVID-19 with piecewise constant dynamical parameters [25] as one possible extension of classical SIR models. Taking age- and sex-structure of a population into account, Wacker and Schlüter proposed an extension of classical SIR models by multiple groups which can be considered as, for example, age groups or locally related groups [26]. Using fractional differential operators, Xu and co-authors analyzed the spread of COVID-19 in 2020 with respect to the United States by means of a fractional-order susceptible-exposed-infected-recovered (SEIR) model [27]. Regarding parameter identification problems in epidemiological models, Marinov and Marinova designed adaptive SIR models to solve these inverse problems [28,29]. In context of inverse problems, Comunian, Gaburro and Giudici concluded that high-quality data are needed to calibrate constant parameters of epidemiological models [30]. Other interesting lines of research concern the consideration of epidemics on networks [31,32,33,34,35] or application of deep-learning approaches to, for example, forecast time-series data of epidemics by these methods [36,37,38].
One important aspect in epidemiological models is preservation of positivity. It plays a huge role in different scientific problems such as wave equations [39], simulations of chemical reactions [7], population dynamics [12,13,14] or other biological systems [3]. However, many classical time discretization schemes suffer from conserving non-negativity for arbitraty time step sizes [40]. Hence, designing non-negativity preserving numerical solution schemes for such problems can be regarded as an art as stated by Mickens [41]. For these reasons, Wacker and Schlüter examined a numerical solution scheme
{Sh,j+1−Sh,jhj+1=−αj+1⋅Sh,j+1⋅Ih,j+1N,Ih,j+1−Ih,jhj+1=αj+1⋅Sh,j+1⋅Ih,j+1N−βj+1⋅Ih,j+1,Rh,j+1−Rh,jhj+1=βj+1⋅Ih,j+1,Sh,1=S0,Ih,1=I0,Rh,1=R0 | (1.2) |
for (1.1) based on implicit Eulerian time discretization in [15] whose algorithm is presented in Algorithm 1.
Algorithm 1 Algorithmic summary for implicit Eulerian finite-difference-method (1.2) for (1.1) |
1: Inputs
2: Time vector (t1,…,tM)T∈RM 3: Time-varying transmission rates (α1,…,αM)T∈RM 4: Time-varying recovery rates (β1,…,βM)T∈RM 5: Initialize S:=0,I:=0,R:=0∈RM 6: Initial conditions S1:=S1, I1:=I1 and R1:=R1 7: Total population size N:=S1+I1+R1 8: Implicit Eulerian Finite-Difference-Method 9: for j∈{1,…,M−1} do 10: Define hj+1:=tj+1−tj as current time step size 11: Define Aj+1:=(1+βj+1⋅hj+1)⋅(αj+1⋅hj+1) 12: Define Bj+1:=(1+βj+1⋅hj+1)⋅N−αj+1⋅hj+1⋅(Sj+Ij)2 13: Compute Ij+1:=−Bj+1Aj+1+√B2j+1A2j+1+N⋅IjAj+1 14: Compute Sj+1:=Sj1+αj+1⋅hj+1⋅Ij+1N 15: Compute Rj+1:=Rj+βj+1⋅hj+1⋅Ij+1 16: end for 17: Outputs 18: Calculated vectors S,I,R∈RM |
Such non-standard finite-difference-methods for time discretization are also important in other application areas such as magnetohydrodynamics [42]. Here, we use f(tj):=fj as an abbreviation for time-continuous functions at time points tj and fh(tj):=fh,j as an abbreviation for time-discrete functions at time time points tj on possibly non-equidistant meshes t1<t2<…<tM−1<tM. Although Algorithm 1 preserves all desired properties of the time-continuous model (1.1) as shown in [15], there are still some disadvantages:
1) From its formulation (1.1), there is to solve a quadratic equation
Aj+1⋅I2j+1+Bj+1⋅Ij+1−N⋅Ij=0 |
where the a-priori knowledge of non-negative solutions from the time-continuous model is used to rule out one of both possible solutions for Algorithm 1 to obtain unique solvability;
2) In each time step, we must recalculate the auxiliary quantities Aj+1 and Bj+1;
3) By both aforementioned arguments, we conclude that Algorithm 1 does not intrinsically calculate biologically correct solutions and, from an algorithm point of view, recalculation of these auxiliary quantities is undesirable.
We want to propose a different discretization by non-standard finite-difference-methods based on non-local approximations. To explain the idea of non-standard finite-difference-methods based on non-local approximations, we consider a time-continuous dynamical system
y′(t)=f(t,y(t)). |
Standard time-discretizations normally approximate the right-hand-side function by function values of y(t) at one previous or the current time point. Non-standard finite-difference-methods based on non-local approximations, on the contrary, approximate f(t,y(t)) by function values from more then solely one previous or the current time point. For further details, we refer to the book [41] by Mickens. For the aforementioned three reasons, we suggest a non-standard finite-difference-method (NSFDM) given by
{Sh,j+1−Sh,jhj+1=−αj+1⋅Sh,j+1⋅Ih,jN,Ih,j+1−Ih,jhj+1=αj+1⋅Sh,j+1⋅Ih,jN−βj+1⋅Ih,j+1,Rh,j+1−Rh,jhj+1=βj+1⋅Ih,j+1,Sh,1=S0,Ih,1=I0,Rh,1=R0 | (1.3) |
in Section 3 as our first main contribution, after summarizing analytical properties of the time-continuous problem formulation (1.1) in Section 2. In Section 3, we additionally prove unconditional preservation of non-negativity and linear convergence of (1.3) towards the time-continuous solution of (1.1). Additionally, (1.3) transfers to a uniquely solvable algorithm which preserves non-negativity intrinsically. Based on (1.3), we design a numerical algorithm to solve the inverse problem of finding the time-varying parameter functions in Section 4 as our next main contribution. Finally, we strengthen our theoretical findings in Section 5 by providing some numerical examples, especially one numerical example where linear convergence is stressed, in contrast to [15], where linear convergence was not shown numerically.
Following [15], we take the following assumptions into consideration with regard to our time-continuous model (1.1):
(1). The total population size N is fixed for all time t≥0, i.e., it holds
N(t)=N=S0+I0+R0 | (2.1) |
for all t≥0 and we assume I0>0 to really get an evolution in time of the infected subgroup of the population;
(2). We divide the total population into three homogeneous subgroups, namely susceptible people (S), infected people (I) and recovered people (R). The equation
N(t)=S(t)+I(t)+R(t) | (2.2) |
is valid for all t≥0;
(3). The time-varying transmission rate α:[0,∞)⟶[0,∞) is at least continuously differentiable once and there exist positive constants αmin>0 and αmax>0 such that
αmin≤α(t)≤αmax | (2.3) |
holds for all t≥0;
(4). The time-varying recovery rate β:[0,∞)⟶[0,∞) is at least continuously differentiable once and there exist positive constants βmin>0 and βmax>0 such that
βmin≤β(t)≤βmax | (2.4) |
holds for all t≥0.
We summarize the following analytical properties of the time-continuous problem formulation (1.1) from article [15]:
Theorem 1. Let all assumptions from (2.1)–(2.4) for our time-continuous model (1.1) be fulfilled. The following statements hold true:
(1). Possible solutions of (1.1) remain non-negative and bounded above, i.e.,
{0≤S(t)≤N;0≤I(t)≤N;0≤R(t)≤N | (2.5) |
hold for all t≥0;
(2). The time-continuous problem formulation (1.1) possesses exactly one globally unique solution for all t≥0;
(3). First, S is monotonically decreasing and there exists S⋆>0 such that limt→∞S(t)=S⋆. Secondly, R is monotonically increasing and there exists R⋆>0 such that limt→∞R(t)=R⋆. Finally, it holds limt→∞I(t)=0.
Proof. (1). This statement's proof can be found in [15, Theorem 4].
(2). Globally unique existence in time is proven in [15, Theorem 6].
(3). Proof of the solution components' monotonicity and its long-time behavior is given in [15, Theorem 8].
We want to note that S⋆ and R⋆ depend on the initial conditions.
Let 0=t1<t2<…<tM=T be an arbitrary and possibly non-equidistant partition of the simulation interval [0,T] with final simulation time T. As abbreviations, we write f(tj):=fj for time-continuous functions and fh(tj):=fh,j for time-discrete functions at time points tj for all j∈{1,2,…,M}. The time-discrete problem formulation of our time-continuous model (1.1) is given by (1.3) for all j∈{1,2,…,M−1}. Here, hj+1:=tj+1−tj denote non-equidistant time step sizes for all j∈{1,2,…,M−1}. Additionally, the supremum norm on [0,T] for time-continuous functions is denoted by
‖f(t)‖∞:=supt∈[0,T]|f(t)| |
while the same supremum norm on [tp,tp+1] reads
‖f(t)‖∞,p+1:=supt∈[tp,tp+1]|f(t)|. |
As our first result, our time-discrete problem formulation (1.3) conserves total population size N for all time points tj∈[0,T] with j∈{1,2,…,M}.
Theorem 2. Possible solutions of (1.3) fulfill
N=Sh,j+Ih,j+Rh,j | (3.1) |
for all j∈{1,2,…,M}.
Proof. This statement can be proven by induction. For j=1, it is clear by our initial conditions that
N=Sh,1+Ih,1+Rh,1 |
holds. Let us assume that
N=Sh,j+Ih,j+Rh,j |
is valid for an arbitrary j∈{1,2,…,M−1}. Adding all three equations of (1.3), we obtain
(Sh,j+1−Sh,j)+(Ih,j+1−Ih,j)+(Rh,j+1−Rh,j)hj+1=0 |
and this implies
Sh,j+1+Ih,j+1+Rh,j+1=Sh,j+Ih,j+Rh,j |
which proves our assertion (3.1).
As our next result, we show that approximate solutions of our time-continuous model (1.1) from our discretization model (1.3) remain non-negative and bounded.
Theorem 3. Possible solutions of (1.3) remain non-negative and bounded for all tj∈[0,T] with arbitrary j∈{1,2,…,M}, i.e., these solultions fulfill
{0≤Sh,j≤N;0≤Ih,j≤N;0≤Rh,j≤N | (3.2) |
for all j∈{1,2,…,M}.
Proof. (1). By the first three equations of (1.3), we obtain
Sh,j+1=Sh,j1+hj+1⋅αj+1⋅Ih,jN | (3.3) |
for all j∈{1,2,…,M−1},
Ih,j+1=Ih,j⋅1+hj+1⋅αj+1Sh,j+1⋅Ih,jN1+hj+1⋅βj+1 | (3.4) |
for all j∈{1,2,…,M−1} and
Rh,j+1=Rh,j+hj+1⋅βj+1⋅Ih,j+1 | (3.5) |
for all j∈{1,2,…,M−1}.
(2). By our assumption of non-negative initial values Sh,1,Ih,1,Rh,1 and by the three equations (3.3)–(3.5), we conclude that
Sh,j≥0;Ih,j≥0;Rh,j≥0 |
hold for all j∈{1,2,…,M}.
(3). By non-negativity and conservation of total population size (3.1) from Theorem 2, it follows that
Sh,j≤N;Ih,j≤N;Rh,j≤N |
hold for all j∈{1,2,…,M}.
(4). Hence, all possible solutions of (1.3) remain non-negative and bounded.
It directly follows that our numerical solution scheme (1.3) is uniquely solvable for all j∈{1,2,…,M}.
Theorem 4. The discretization scheme (1.3), based on the explicit solutions (3.3)–(3.5), is well-defined and uniquely solvable for all j∈{1,2,…,M}.
Proof. (3.3)–(3.5) prove that our numerical solution scheme (1.3) is well-defined. Unique solvability is a direct consequence of these equations as well. Hence, our assertion is proven.
Remark 1. We briefly want to note why we call solutions (3.3)–(3.5) explicit. It obviously holds
Sh,j+1=Sh,j1+hj+1⋅αj+1⋅Ih,jN=N⋅Sh,jN+hj+1⋅αj+1⋅Ih,j | (3.6) |
for all j∈{1,2,…,M} explicitly. Plugging (3.6) into (3.4) yields
Ih,j+1=Ih,j⋅(1+hj+1⋅αj+1⋅Sh,j+1⋅Ih,jN)1+hj+1⋅βj+1=Ih,j⋅(1+hj+1⋅αj+1⋅Sh,j⋅Ih,jN+hj+1⋅αj+1⋅Ih,j)1+hj+1⋅βj+1 | (3.7) |
for all j∈{1,2,…,M} as an explicit formula. Finally, by using (3.7), we get
Rh,j+1=Rh,j+hj+1⋅βj+1⋅Ih,j+1=Rh,j+hj+1⋅βj+1⋅Ih,j⋅(1+hj+1⋅αj+1⋅Sh,j⋅Ih,jN+hj+1⋅αj+1⋅Ih,j)1+hj+1⋅βj+1 | (3.8) |
for all j∈{1,2,…,M}. For these reasons, we talk about an explicit numerical solution scheme.
We show that monotonicity properties and long-time behavior of the time-continuous problem formulation (1.1) translates to our time discretization (1.3).
Theorem 5. 1). The sequence {Sh,j}Mj=1 is monotonically decreasing and there exists a non-negative constant S⋆h such that limj→∞Sh,j=S⋆h holds.
2). The sequence {Rh,j}Mj=1 is monotonically increasing and there exists a non-negative constant R⋆h such that limj→∞Rh,j=R⋆h holds.
3). It holds limj→∞Ih,j=0 if also T→∞ and hj+1≠0 are assumed.
Proof. (1). Since {Sh,j}Mj=1 is non-negative and bounded above by the total population size N, the inequality
Sh,j+1=Sh,j1+hj+1⋅αj+1⋅Ih,jN≤Sh,j |
from (3.6) implies that this sequence is also monotonically decreasing which proves our first claim.
(2). Since {Rh,j}Mj=1 is non-negative and bounded above by the total population size N, the inequality
Rh,j+1=Rh,j+hj+1⋅βj+1⋅Ij+1≥Rh,j |
from (3.5) yields that this sequence is monotonically increasing which shows our second assertion.
(3). We conclude
0≤Ij+1=Rh,j+1−Rh,jhj+1⋅βj+1≤Rh,j+1−Rh,jhj+1⋅βmin |
from (3.5) and this shows that our third claim is also true.
Remark 2. Both Theorems 1 and 5 imply that our time-continuous problem formulation (1.1) and our time-discrete problem formulation (1.3) converge towards the disease-free equilibrium state in the long run.
We start this subsection by listing all our assumptions for our main result regarding linear convergence of the time-discrete solution towards the time-continuous one:
(A11) We consider the time interval [0,T] where
t1=0<t2<…<tM−1<tM:=T |
holds;
(A12) The initial conditions of the time-continuous and of the time-discrete problems coincide;
(A13) The time-continuous solution functions S,I,R:[0,T]⟶R should be twice continuously differentiable;
(A14) The time-varying transmission rate α:[0,T]⟶[0,∞) and the time-varying recovery rate β:[0,T]⟶[0,∞) should be once continuously differentiable;
(A15) The time-varying transmission rate and the time-varying recovery rate should be bounded, i.e., there are non-negative constants αmin,αmax,βmin,βmax such that
0≤αmin≤α(t)≤αmax |
and
0≤βmin≤β(t)≤βmax |
hold for all t∈[0,T];
(A16) Set h:=maxp∈Nhp.
Theorem 6. Let all assumptions (A11)–(A16) be fulfilled. Then the time-discrete solution convergences linearly towards the time-continuous solution if h→0.
Proof. We shortly explain our proof's strategy because it is relatively technical. If we first assume that our time-discrete solution coincides with the time-continuous solution at a certain time point tp∈[0,T] for an arbitrary p∈{1,2,…,M−1}, then we investigate the local error propagation to the next time point tp+1. Secondly, we consider error propagation in time and finally, we examine accumulation of these errors. We denote time-discrete solutions, for example, by Sh,p and time-continuous solutions by S(tp) at the same time point tp.
1) For investigation of local errors, we assume that
(tp,Sh,p)=(tp,S(tp));(tp,Ih,p)=(tp,I(tp))and(tp,Rh,p)=(tp,R(tp)) |
hold for an arbitrary p∈{1,2,…,M−1} on the time-interval [tp,tp+1]. Since we only take one time step into account, we denote corresponding time-discrete solutions by ~Sh,p+1, ~Ih,p+1 and ~Rh,p+1.
1.1) We have by (3.6) from Remark 1 that
~Sh,p+1=S(tp)(1+αp+1⋅hp+1⋅I(tp)N)=S(tp)−S(tp)⋅αp+1⋅hp+1⋅I(tp)N(1+αp+1⋅hp+1⋅I(tp)N) | (3.9) |
holds. This implies
|S(tp+1)−~Sh,p+1|(3.9)=|S(tp+1)−S(tp)+S(tp)⋅αp+1⋅hp+1⋅I(tp)N(1+αp+1⋅hp+1⋅I(tp)N)|≤|tp+1∫tpS′(τ)dτ−hp+1⋅S′(tp)|⏟=:IS,1+|hp+1⋅S′(tp)+S(tp)⋅αp+1⋅hp+1⋅I(tp)N(1+αp+1⋅hp+1⋅I(tp)N)|⏟=:IIS,1. | (3.10) |
By the mean value Theorem, there exists ξS,1∈(tp,tp+1) such that we obtain
|S′′(ξS,1)|=|S′(τ)−S′(tp)τ−tp|≤‖S′′(t)‖∞. | (3.11) |
This yields
IS,1=|tp+1∫tp{S′(τ)−S′(tp)}dτ|=|tp+1∫tp(S′(τ)−S′(tp))(τ−tp)⋅(τ−tp)dτ|(3.11)≤12⋅h2p+1⋅‖S′′(t)‖∞. | (3.12) |
For IIS,1, we notice that
IIS,1=|hp+1⋅S′(tp)+S(tp)⋅αp+1⋅hp+1⋅I(tp)N(1+αp+1⋅hp+1⋅I(tp)N)|(1.1)=|−αp⋅hp+1⋅S(tp)⋅I(tp)N+S(tp)⋅αp+1⋅hp+1⋅I(tp)N(1+αp+1⋅hp+1⋅I(tp)N)|=|hp+1⋅S(tp)N|⋅|−αp⋅I(tp)+αp+1⋅I(tp)(1+αp+1⋅hp+1⋅I(tp)N)|≤hp+1⋅|−αp⋅I(tp)⋅(1+αp+1⋅hp+1⋅I(tp)N)+αp+1⋅I(tp)(1+αp+1⋅hp+1⋅I(tp)N)|≤hp+1⋅|(αp+1−αp)⋅I(tp)(1+αp+1⋅hp+1⋅I(tp)N)|+hp+1⋅|αp⋅αp+1⋅hp+1⋅I(tp)⋅I(tp)N(1+αp+1⋅hp+1⋅I(tp)N)|≤hp+1⋅N⋅|αp+1−αp|⏟=:IIIS,1+h2p+1⋅α2max⋅N |
and hence,
IIS,I≤hp+1⋅N⋅|αp+1−αp|⏟=:IIIS,1+h2p+1⋅α2max⋅N | (3.13) |
holds due to boundedness of time-continuous solutions. Using again the mean value theorem for IIIS,1, we see that there is ξα,1∈(tp,tp+1) such that we obtain
|α′(ξα,1)|=|αp+1−αptp+1−tp|⋅hp+1≤hp+1⋅‖α′(t)‖∞. | (3.14) |
Hence, we conclude
IIIS,1=|(αp+1−αp)|(3.14)≤hp+1⋅‖α′(t)‖∞ | (3.15) |
for IIIS,1. Plugging (3.15) into (3.13) yields
IIS,1≤hp+1⋅N⋅IIIS,1+h2p+1⋅α2max⋅N(3.15)≤h2p+1⋅N⋅‖α′(t)‖∞+h2p+1⋅α2max⋅N=h2p+1⋅(‖α′(t)‖∞⋅N+α2max⋅N). | (3.16) |
Thus, we finally obtain
|S(tp+1)−~Sh,p+1|(3.10)=IS,1+IIS,1≤12⋅h2p+1⋅‖S′′(t)‖∞+h2p+1⋅(‖α′(t)‖∞⋅N+α2max⋅N)=(12⋅‖S′′(t)‖∞+‖α′(t)‖∞⋅N+α2max⋅N)⏟=:Cloc,S⋅h2p+1. | (3.17) |
1.2) By (3.7) from Remark 1, we get
~Ih,p+1=I(tp)⋅(1+αp+1⋅hp+1⋅~Sh,p+1N)(1+βp+1⋅hp+1). | (3.18) |
Hence, it follows
|I(tp+1)−~Ih,p+1|(3.18)=|I(tp+1)−I(tp)⋅(1+αp+1⋅hp+1⋅~Sh,p+1N)(1+βp+1⋅hp+1)|=|{I(tp+1)−I(tp)}+{I(tp)−I(tp)⋅(1+αp+1⋅hp+1⋅~Sh,p+1N)(1+βp+1⋅hp+1)}|=|tp+1∫tpI′(τ)dτ+I(tp)⋅(βp+1⋅hp+1−αp+1⋅hp+1⋅~Sh,p+1N)(1+βp+1⋅hp+1)|≤|tp+1∫tpI′(τ)dτ−hp+1⋅I′(tp+1)|⏟=:II,1+|hp+1⋅I′(tp+1)+I(tp)⋅(βp+1⋅hp+1−αp+1⋅hp+1⋅~Sh,p+1N)(1+βp+1⋅hp+1)|⏟=:III,1. | (3.19) |
By the mean value theorem, there exists ξI,1∈(tp,tp+1) such that
I′′(ξI,1)=I′(τ)−I′(tp+1)τ−tp+1 | (3.20) |
is valid. This implies
II,1=|tp+1∫tp(I′(τ)−I′(tp+1))(τ−tp+1)⋅(τ−tp+1)dτ|(3.20)=|tp+1∫tpI′′(ξI,1)⋅(τ−tp+1)dτ|≤12⋅‖I′′(t)‖∞⋅h2p+1. | (3.21) |
For III,1, we observe
III,1=|hp+1⋅I′(tp+1)+I(tp)⋅(βp+1⋅hp+1−αp+1⋅hp+1⋅~Sh,p+1N)(1+βp+1⋅hp+1)|=hp+1⋅|αp+1⋅I(tp+1)⋅S(tp+1)N−βp+1⋅I(tp+1)+I(tp)⋅(βp+1⋅hp+1−αp+1⋅hp+1⋅~Sh,p+1N)(1+βp+1⋅hp+1)|=hp+11+βp+1⋅hp+1⋅|αp+1I(tp+1)⋅S(tp+1)N⋅(1+βp+1⋅hp+1)−βp+1⋅I(tp+1)⋅(1+βp+1⋅hp+1)+βp+1⋅I(tp)−αp+1⋅I(tp)⋅~Sh,p+1N| | (3.22) |
and can conclude
III,1(3.22)≤hp+1⋅|αp+1N⋅(I(tp+1)⋅S(tp+1)−I(tp)⋅~Sh,p+1)|⏟=:IIII,1+hp+1⋅|αp+1⋅βp+1⋅hp+1⋅I(tp+1)⋅S(tp+1)N|⏟=:IVI,1+hp+1⋅|βp+1⋅(I(tp)−I(tp+1))|⏟=:VI,1+hp+1⋅|β2p+1⋅I(tp+1)⋅hp+1|⏟=:VII,1. | (3.23) |
By boundedness and the mean value theorem, we receive the following estimates
VII,1≤hp+1⋅N⋅β2max,VI,1≤hp+1⋅βmax⋅‖I′(t)‖∞,IVI,1≤hp+1⋅N⋅αmax⋅βmax | (3.24) |
for VII,1, VI,1 and IVI,1. By boundedness, the triangle inequality and the mean value theorem, we can deduce
IIII,1=|αp+1N|⋅|{I(tp+1)⋅S(tp+1)−I(tp)⋅S(tp+1)}+{I(tp)⋅S(tp+1)−I(tp)⋅~Sh,p+1}|≤αmaxN⋅|I(tp+1)−I(tp)|⋅|S(tp+1)|+αmaxN⋅|S(tp+1)−~Sh,p+1|⋅|I(tp)|(3.17)≤αmax⋅hp+1⋅‖I′(t)‖∞+αmax⋅h2p+1⋅Cloc,S | (3.25) |
and consequently, this implies
III,1≤hp+1⋅{IIII,1+IVI,1+VI,1+VII,1}≤hp+1⋅{αmax⋅hp+1⋅‖I′(t)‖∞+αmax⋅h2p+1⋅Cloc,S+αmax⋅βmax⋅hp+1⋅N+βmax⋅hp+1⋅‖I′(t)‖∞+β2max⋅hp+1⋅N}. | (3.26) |
Hence, we conclude that
|I(tp+1)−~Ih,p+1|(3.19)≤II,1+III,1(3.26)≤12⋅h2p+1⋅‖I′′(t)‖∞+h2p+1⋅{αmax⋅‖I′(t)‖∞+αmax⋅hp+1⋅Cloc,S+αmax⋅βmax⋅N+βmax⋅‖I′(t)‖∞+β2max⋅N} | (3.27) |
holds. By defining
Cloc,I:={12⋅‖I′′(t)‖∞+αmax⋅‖I′(t)‖∞+αmax⋅hp+1⋅Cloc,S+αmax⋅βmax⋅N+βmax⋅‖I′(t)‖∞+β2max⋅N}, | (3.28) |
we finally deduce
|I(tp+1)−~Ih,p+1|(3.28)≤h2p+1⋅Cloc,I | (3.29) |
from (3.27).
1.3) By (3.8) from Remark 1, we get
~Rh,p+1=R(tp)+hp+1⋅βp+1⋅~Ih,p+1. | (3.30) |
Hence, it follows
|R(tp+1)−~Rh,p+1|(3.30)=|R(tp+1)−R(tp)−hp+1⋅βp+1⋅~Ih,p+1|≤|tp+1∫tpR′(τ)dτ−hp+1⋅βp+1⋅I(tp+1)|+hp+1⋅βmax⋅|I(tp+1)−~Ih,p+1|. | (3.31) |
By using similar argument as in previous steps with the help of the mean value theorem and (3.29), we can deduce
|R(tp+1)−~Rh,p+1|≤12⋅h2p+1⋅‖R′′(t)‖∞+h3p+1⋅βmax⋅Cloc,I=h2p+1⋅{12⋅‖R′′(t)‖∞+βmax⋅Cloc,I⋅hp+1}⏟=:Cloc,R. | (3.32) |
1.4) By setting
Cloc:=max{Cloc,S,Cloc,I,Cloc,R}, | (3.33) |
we conclude
Ap+1:=max{|S(tp+1)−~Sh,p+1|,|I(tp+1)−~Ih,p+1|,|R(tp+1)−~Rh,p+1|}≤Cloc⋅h2p+1. | (3.34) |
2) Normally, the points (tp,Sp), (tp,Ip) and (tp,Rp) do not exactly lie on the graph of the time-continuous solution. We then must investigate how procedural errors such as |Sh,p−S(tp)|, |Ih,p−I(tp)| or |Rh,p−R(tp)| propagate to the (p+1)-th time step.
2.1) We observe that
|Sh,p+1−~Sh,p+1|=|(Sh,p−Sh,p⋅αp+1⋅hp+1⋅Ih,pN1+αp+1⋅hp+1⋅Ih,pN)−(S(tp)−S(tp)⋅αp+1⋅hp+1⋅I(tp)N1+αp+1⋅hp+1⋅I(tp)N)|≤|Sh,p−S(tp)|+αmax⋅hp+1N⋅IS,2 | (3.35) |
holds where IS,2 is given and estimated by
IS,2=|Sh,p⋅Ih,p⋅(1+αp+1⋅hp+1⋅I(tp)N)−S(tp)⋅I(tp)⋅(1+αp+1⋅hp+1⋅Ih,pN)(1+αp+1⋅hp+1⋅Ih,pN)⋅(1+αp+1⋅hp+1⋅I(tp)N)|≤|Sh,p⋅Ih,p−S(tp)⋅I(tp)|⏟=:IIS,2+αmax⋅hp+1N⋅|Sh,p⋅Ih,p⋅I(tp)−S(tp)⋅I(tp)⋅Ih,p|⏟=:IIIS,2. | (3.36) |
We have
IIS,2=|Sh,p⋅Ih,p−S(tp)⋅I(tp)|≤|Sh,p⋅(Ih,p−I(tp))|+|I(tp)⋅(Sh,p−S(tp))|≤N⋅|Ih,p−I(tp)|+N⋅|Sh,p−S(tp)| | (3.37) |
and
IIIS,2=|Sh,p⋅Ih,p⋅I(tp)−S(tp)⋅I(tp)⋅Ih,p|=|I(tp)⋅(Sh,p⋅Ih,p−S(tp)⋅Ih,p)|≤N⋅|Ih,p⋅(Sh,p−S(tp))|≤N2⋅|Sh,p−S(tp)|. | (3.38) |
Plugging (3.37) and (3.38) into (3.36) yields
IS,2≤IIS,2+αmax⋅hp+1N⋅IIIS,2≤N⋅|Ih,p−I(tp)|+N⋅|Sh,p−S(tp)|+αmax⋅hp+1⋅N⋅|Sh,p−S(tp)|=N⋅|Ih,p−I(tp)|+N⋅(1+αmax⋅hp+1)⋅|Sh,p−S(tp)|. | (3.39) |
Using (3.35), we obtain
|Sh,p+1−~Sh,p+1|≤|Sh,p−S(tp)|+αmax⋅hp+1N⋅IS,2≤(1+αmax⋅hp+1+(αmax⋅hp+1)2)⋅|Sh,p−S(tp)|+αmax⋅hp+1⋅|Ih,p−I(tp)|. | (3.40) |
2.2) We see that
|Ih,p+1−~Ih,p+1|=|(Ih,p−Ih,p⋅(βp+1⋅hp+1−αp+1⋅hp+1⋅Sh,p+1N)1+βp+1⋅hp+1)−(I(tp)−I(tp)⋅(βp+1⋅hp+1−αp+1⋅hp+1⋅~Sh,p+1N)1+βp+1⋅hp+1)| | (3.41) |
holds and get
|Ih,p+1−~Ih,p+1|≤|Ih,p−I(tp)|+βmax⋅hp+1⋅|Ih,p−I(tp)|+αmax⋅hp+1N⋅|Ih,p⋅Sh,p+1−I(tp)⋅~Sh,p+1|≤|Ih,p−I(tp)|+βmax⋅hp+1⋅|Ih,p−I(tp)|+αmax⋅hp+1N⋅|Ih,p⋅Sh,p+1−I(tp)⋅Sh,p+1|+αmax⋅hp+1N⋅|I(tp)⋅Sh,p+1−I(tp)⋅~Sh,p+1|≤(1+hp+1⋅(αmax+βmax))⋅|Ih,p−I(tp)|+αmax⋅hp+1⋅|Sh,p+1−~Sh,p+1|. | (3.42) |
Plugging (3.40) into (3.42) yields
|Ih,p+1−~Ih,p+1|≤(1+hp+1⋅(αmax+βmax))⋅|Ih,p−I(tp)|+αmax⋅hp+1⋅{(1+αmax⋅hp+1+(αmax⋅hp+1)2)⋅|Sh,p−S(tp)|+αmax⋅hp+1⋅|Ih,p−I(tp)|}=(1+hp+1⋅(αmax+βmax)+α2max⋅h2p+1)⋅|Ih,p−I(tp)|+(αmax⋅hp+1+(αmax⋅hp+1)2+(αmax⋅hp+1)3)⋅|Sh,p−S(tp)|. | (3.43) |
2.3) We obtain
|Rh,p+1−~Rh,p+1|≤|Rh,p−R(tp)|+βmax⋅hp+1⋅|Ih,p+1−~Ih,p+1|(3.43)≤|Rh,p−R(tp)|+(βmax⋅hp+1+βmax⋅h2p+1⋅(αmax+βmax)+βmax⋅α2max⋅h3p+1)⋅|Ih,p−I(tp)|+(βmax⋅αmax⋅h2p+1+βmax⋅α2max⋅h3p+1+βmax⋅α3max⋅h4p+1)⋅|Sh,p−S(tp)| | (3.44) |
2.4) For our summary, we define
Bp+1:=max{|Sh,p+1−~Sh,p+1|,|Ih,p+1−~Ih,p+1|,|Rh,p+1−~Rh,p+1|} | (3.45) |
for all p∈{0,1,…,M−1} and
Dp:=max{|Sh,p−S(tp)|,|Ih,p−I(tp)|,|Rh,p−R(tp)|} | (3.46) |
for all p∈{1,…,M}. By (3.40), (3.43) and (3.44), we define the constants
CS,prop:=1+hp+1⋅2⋅αmax+h2p+1⋅α2max=1+hp+1⋅(2⋅αmax+hp+1⋅α2max)⏟=:~CS,prop;CI,prop:=1+hp+1⋅(2⋅(αmax+βmax)+2⋅α2max⋅hp+1+α3max⋅h2p+1)⏟=:~CI,prop;CR,prop:=1+hp+1⋅(βmax+2⋅βmax⋅(αmax+βmax)⋅hp+1+2⋅βmax⋅α2max⋅h2p+1+βmax⋅α3max⋅h3p+1)=:1+hp+1⋅~CR,prop |
and
Cprop:=max{~CS,prop,~CI,prop,~CR,prop}. | (3.47) |
Hence, it holds
Bp+1≤(1+hp+1⋅Cprop)⋅Dp. | (3.48) |
3) We obtain
Dp+1≤Bp+1+Ap+1≤(1+hp+1⋅Cprop)⋅Dp+Cloc⋅h2p+1≤(1+hp+1⋅Cprop)⋅((1+hp⋅Cprop)⋅Dp−1+Cloc⋅h2p)+Cloc⋅h2p+1 | (3.49) |
and inductively
Dp+1≤exp(Cprop⋅T)⋅D1+exp(Cprop⋅T)⋅Cloc⋅h⋅T | (3.50) |
which proves linear convergence of the time-discrete solution towards the time-continuous solution.
Our algorithm for our proposed non-standard finite-difference-method for (1.1) is portrayed in Algorithm 2 below.
Algorithm 2 Algorithmic summary for our proposed non-standard finite-difference-method (3.2) for (1.1) |
1: Inputs
2: Time vector (t1,…,tM)T∈RM 3: Time-varying transmission rates (α1,…,αM)T∈RM 4: Time-varying recovery rates (β1,…,βM)T∈RM 5: Initialize S:=0,I:=0,R:=0∈RM 6: Initial conditions S1:=S1, I1:=I1 and R1:=R1 7: Total population size N:=S1+I1+R1 8: Non-Standard Finite-Difference-Method 9: for j∈{1,…,M−1} do 10: hj+1:=tj+1−tj 11: Sj+1:=Sj1+hj+1⋅αj+1⋅IjN from (3.6) 12: Ij+1:=Ij+hj+1⋅αj+1⋅Sj+1⋅IjN1+hj+1⋅βj+1 from (3.7) 13: Rj+1:=Rj+hj+1⋅βj+1⋅Ij+1 from (3.8) 14: end for 15: Outputs 16: Calculated vectors S,I,R∈RM |
We consider the parameter identification problem (1.3) for all j∈{1,2,…,M−1} where we take the following assumptions into account:
(A21) The sequences {Sh,p}Mp=1, {Ih,p}Mp=1 and {Rh,p}Mp=1 are given;
(A22) It holds 0≤Sh,p≤N for all p∈{1,…,M} and the sequence {Sh,p}Mp=1 is monotonically decreasing;
(A23) It holds 0≤Ih,p≤N for all p∈{1,…,M};
(A24) It holds 0≤Rh,p≤N for all p∈{1,…,M} and the sequence {Rh,p}Mp=1 is monotonically increasing;
(A25) We assume Sh,p≠0, Ih,p≠0 and Rh,p≠0 for all p∈{1,…,M}.
Our idea now is to reduce (1.3) to
{Sh,p+1−Sh,php+1=−αp+1⋅Sh,p+1⋅Ih,pN;Rh,p+1−Rh,php+1=βp+1⋅Ih,p+1 | (4.1) |
for all p∈{1,…,M−1}. Hence, we obtain the sought time-varying transmission rate and the sought time-varying recovery rate by
{αp+1=NSh,p+1⋅Ih,p⋅(Sh,p−Sh,p+1)hp+1;βp+1=1Ih,p+1⋅(Rh,p+1−Rh,p)hp+1 | (4.2) |
for all p∈{1,…,M−1}. Finally, we have explicit formulas for the time-varying transmission and recovery rates at every time point which we can evaluate by the known data from the numbers of susceptible, infected and recovered people. We have the following theorem.
Theorem 7. Let all assumptions (A21)–(A25) be fulfilled. (4.1) for parameter identification is uniquely solvable for all p∈{1,…,M−1} and it holds αp+1≥0 and βp+1≥0 for all p∈{1,…,M−1}.
Proof. 1) Unique solvability is a direct consequence of (4.2).
2) {Sh,p}Mp=1 is monotonically decreasing. Hence, Sh,p−Sh,p+1≥0 for all p∈{1,…,M−1} and we have
αp+1=NSh,p+1⋅Ih,p⋅(Sh,p−Sh,p+1)hp+1≥0. |
3) {Rh,p}Mp=1 is monotonically increasing. Hence, Rh,p+1−Rh,p≥0 for all p∈{1,…,M−1} and we have
βp+1=1Ih,p+1⋅(Rh,p+1−Rh,p)hp+1≥0. |
This finishes our claim's proof.
Our algorithmic summary for our proposed parameter identification problem (4.2) reads:
Algorithm 3 Algorithmic summary for our proposed parameter identification algorithm for (4.2) |
1: Inputs
2: Time vector (t1,…,tM)T∈RM 3: Sequence of susceptible population {Sh,p}Mp=1 saved in one vector S∈RM 4: Sequence of infected population {Ih,p}Mp=1 saved in one vector I∈RM 5: Sequence of recovered population {Rh,p}Mp=1 saved in one vector R∈RM 6: Constant total population size N:=S1+I1+R1 7: Initialize vectors α:=0∈RM and β:=0∈RM 8: Parameter Identification Algorithm 9: for j∈{1,…,M−1} do 10: hj+1:=tj+1−tj 11: αj+1:=NSh,j+1⋅Ih,j⋅(Sh,j−Sh,j+1)hj+1 12: βj+1=1Ih,j+1⋅(Rh,j+1−Rh,j)hj+1 13: end for 14: Outputs 15: Calculated vectors α,β∈RM |
In this section, we want to strengthen our theoretical findings by two numerical examples. We apply GNU OCTAVE for our implementation [43].
The initial conditions for our first example read S(0)=5000, I(0)=5000, R(0)=0, N=10,000, α=0.5, β=0.1 and our simulation interval is given by [0,T] with T=70. We compare three numerical algorithms, namely our proposed non-standard finite-difference-method (NSFDM), the explicit Eulerian scheme (EE) and an explicit Runge-Kutta scheme with two stages (RK2) and use equidistant time step sizes h for our computations.
As it can be clearly seen in Figures 1 and 2, only our novel non-standard finite-difference-method is unconditionally non-negativity-preserving and stable. Both explicit standard time-discretizations suffer from time-step restrictions to obtain non-negativity preservation. In Figure 3, we observe that all three numerical algorithms possess same properties for smaller time step sizes h. On the positive side for all algorithms, they conserve total population size N as depicted in Figure 4.
Finally, underlining our convergence result from Theorem 6, we compare results of our non-standard finite-difference-method (NSFDM) for different time step sizes h with a fine scale Runge-Kutta scheme of order 4 (RK4) at T=70 where a fine time step size of 0.0001 is applied for (RK4). The errors
● Error for S(t): |SNSFDM,M−SRK4,M|;
● Error for I(t): |INSFDM,M−IRK4,M|;
● Error for R(t): |RNSFDM,M−RRK4,M|
are used for this purpose. It can be clearly seen from Table 1 that the theoretical convergence order is achieved by our non-standard finite-difference-method.
h | 0.03125 | 0.0625 | 0.125 | 0.25 | 0.5 | 1.0 |
Error for S(t) | 0.27673 | 0.55049 | 1.0889 | 2.1307 | 4.0817 | 7.5129 |
Error for I(t) | 0.14805 | 0.2971 | 0.59806 | 1.2115 | 2.4843 | 5.2145 |
Error for R(t) | 0.41436 | 0.83717 | 1.6766 | 3.3317 | 6.5556 | 12.717 |
We consider real-world data from Spain [44]. N=47,370,000 is the total population size of Spain. The cumulative cases of infected people are modified such
S(t)=N−I(t)−R(t);I(t)=Nconfirmed(t)−Ndead(t)−Nrecovered(t);R(t)=Ndead(t)+Nrecovered(t) |
hold from the given data. We take α(t)=0.52⋅e−0.03⋅t and β(t)=0.045. For our simulation, we additionally assume h=0.75 as an equidistant time step size.
Our estimated transmission rate α(t) and constant recovery rate β(t) are user-chosen and adapted to the data in Figures 6 and 7, while the modified data for infected and recovered people are depicted in Figure 5. We solely consider the beginning of the pandemic because, later, vaccines have been available and reinfections have occurred after virus modifications. In Figure 8, we see results of our numerical simulations with h=0.75. Results are relatively sensitive with respect to time step sizes and estimates rates, e.g., time-varying transmission rates. Hence, we can only conclude that good data acquisition and different predictions in different settings are crucial for future pandemics to make predictions for courses of such future epidemics.
In this article, we first summarized some analytical properties of a time-continuous epidemiological model (1.1) as given in [15]. Afterwards, we designed a new non-standard finite-difference-method (1.3) as a time-discrete version of the continuous model. In Section 3, we contributed our main results on our proposed non-standard finite-difference-method (1.3). Namely, we investigated unconditionally non-negativity preservation, conservation of total population size, convergence towards a disease-free equilibrium point and linear convergence by letting the time size h→0. Conclusively for our theoretical findings, we developped a simple parameter identification algorithm in Section (4). In Section 5, we gave two numerical examples to underline our theoretical findings. Comparing Algorithms 1 and 2, we conclude that the latter one is definitely easier to implement and it can serve as a starting point for constructing higher order time discretization schemes for (1.1). Regarding practical applications, our proposed time-discretization by a non-standard finite-difference-method preserves all desired properties from the time-continuous model which are expected from our real-word intuition. Additionally, the simple parameter identification algorithm produce fast graphs from given data where trends in time-varying transmission rates such as contact restrictions might be seen from.
There might be different directions to extend our work:
● It seems straightforward to integrate further compartments into our model such that we can additionally investigate influences such as, for example, vaccines or the simple extension to the SEIR model. For example, we readily can construct
{Sh,j+1−Sh,jhj+1=−αj+1⋅Sh,j+1⋅Eh,jN,Eh,j+1−Eh,jhj+1=αj+1⋅Sh,j+1⋅Eh,jN−γj+1⋅Eh,j+1,Ih,j+1−Ih,jhj+1=γj+1⋅Eh,j+1−βj+1⋅Ih,j+1,Rh,j+1−Rh,jhj+1=βj+1⋅Ih,j+1,Sh,1=S0,Eh,1=E0,Ih,1=I0,Rh,1=R0 | (6.1) |
as one variant of a non-standard finite-difference-method for the SEIR model and it follows directly that this scheme conserves total population size N and can be written in an explicit fashion as our proposed non-standard finite-difference-method for our considered SIR model by
{Sh,j+1=Sh,j1+hj+1⋅αj+1⋅Eh,jN,Eh,j+1=Eh,j+hj+1⋅αj+1⋅Sh,j+1⋅Eh,jN1+hj+1⋅γj+1,Ih,j+1=Ih,j+hj+1⋅γj+1⋅Eh,j+11+hj+1⋅βj+1,Rh,j+1=Rh,j+hj+1⋅βj+1⋅Ih,j+1 | (6.2) |
to see that it also preserves non-negativity assuming positive initial conditions. One disadvantage is that there are different possibilities to construct such non-standard finite-difference-methods by non-local approximations and a suitable choice depends on the time-continuous properties which should be transferred to the time-discrete case. However, as seen by this example, it is flexible to construct different possibilities of algorithms and they are easy to implement;
● It might be of interest to construct higher order non-standard finite-difference-methods from our first-order non-standard finite-difference-method in order to get higher accuracy although, especially, our second example demonstrated that predictions of such models need to be considered with care. One straightforward way to improve accuracy, and hence, the convergence order of numerical solution schemes is application of Richardson's extrapolation [45];
● One can use our discretization model for the discretization of optimal control problems for minimizing certain epidemic aspect such as eradication time [46].
Our code and used data from [44] can be found under https://github.com/bewa87/2023-SIR-FirstOrder-NSFDM.
We thank the editor and both anonymous reviewers for their valuable comments which definitely helped us to improve quality and readability of our manuscript.
The authors declare there is no conflict of interest.
[1] | Kanis JA, Kanis JA (1994) Assessment of fracture risk and its application to screening for postmenopausal osteoporosis: synopsis of a WHO report. Osteoporosis Int 4: 368–381. |
[2] | Beck GR, Ha SW, Camalier CE, et al. (2012) Bioactive silica-based nanoparticles stimulate bone-forming osteoblasts, suppress bone-resorbing osteoclasts, and enhance bone mineral density in vivo. Nanomed Nanotechnol Biol Med 8: 793–803. |
[3] | Khajuria DK, Razdan R, Mahapatra DR (2011) Drugs for the management of osteoporosis: a review. Rev Bras Reumatol 51: 365–382. |
[4] | Arcos D, Boccaccini AR, Bohner M, et al. (2014) The relevance of biomaterials to the prevention and treatment of osteoporosis. Acta Biomater 10: 1793–1805. |
[5] | Wei D, Jung J, Yang H, et al. (2016) Nanotechnology treatment options for osteoporosis and its corresponding consequences. Curr Osteoporosis Rep 2016: 1–9. |
[6] | Riggs BL, Melton LJ (2015) The worldwide problem of osteoporosis: insights afforded by epidemiology. Bone 17: 505S–511S. |
[7] | Mackey PA, Whitaker MD (2015) Osteoporosis: a therapeutic update. J Nurse Pract 11: 1011–1017. |
[8] | Hernlund E, Svedbom A, Ivergård M, et al. (2013) Osteoporosis in the European Union: medical management, epidemiology and economic burden: a report prepared in collaboration with the international osteoporosis foundation (IOF) and the european federation of pharmaceutical industry associations (EFPIA). Arch Osteoporosis 8: 1–115. |
[9] | Kanis JA, McCloskey EV, Johansson H, et al. (2013) European guidance for the diagnosis and management of osteoporosis in postmenopausal women. Osteoporosis Int 24: 23–57. |
[10] | Kim T, Singh RK, Sil M, et al. (2016) Acta biomaterialia inhibition of osteoclastogenesis through siRNA delivery with tunable mesoporous bioactive nanocarriers. Acta Biomater 29: 352–364. |
[11] | Hodsman AB, Bauer DC, Dempster DW, et al. (2005) Parathyroid hormone and teriparatide for the treatment of osteoporosis: a review of the evidence and suggested guidelines for its use. Endocr Rev 26: 688–703. |
[12] | Weitzmann MN, Ha SW, Vikulina T, et al. (2015) Bioactive silica nanoparticles reverse age-associated bone loss in mice. Nanomed Nanotechnol Biol Med 11: 959–967. |
[13] | Khajuria DK, Disha C, Vasireddi R, et al. (2016) Risedronate/zinc-hydroxyapatite based nanomedicine for osteoporosis. Mat Sci Eng C Mater 63: 78–87. |
[14] | Giger EV, Castagner B, Leroux JC (2013) Biomedical applications of bisphosphonates. J Control Release 167: 175–188. |
[15] | Iñiguez ANM, Clarke BL (2015) Bone biology, signaling pathways, and therapeutic targets for osteoporosis. Maturitas 82: 245–255. |
[16] | Kanakaris NK, Petsatodis G, Tagil M, et al. (2009) Is there a role for bone morphogenetic proteins in osteoporotic fractures? Injury 40: S21–S26. |
[17] | Lee D, Heo DN, Kim H, et al. (2016) Inhibition of osteoclast differentiation and bone resorption by bisphosphonate- conjugated gold nanoparticles. Sci Rep 6: 27336–27346. |
[18] | Eriksen EF, Díez PA, Boonen S (2014) Update on long-term treatment with bisphosphonates for postmenopausal osteoporosis: a systematic review. Bone 58: 126–135. |
[19] | Black DM, Schwartz AV, Ensrud KE, et al. (2006) Effects of continuing or stopping alendronate after 5 years of treatment: the fracture intervention trial long-term extension (FLEX): a randomized trial. JAMA 296: 2927–2938. |
[20] | Dougall WC, Glaccum M, Charrier K, et al. (1999) RANK is essential for osteoclast and lymph node development. Gene Dev 13: 2412–2424. |
[21] | Trejo CG, Lozano D, Manzano M, et al. (2010) The osteoinductive properties of mesoporous silicate coated with osteostatin in a rabbit femur cavity defect model. Biomaterials 31: 8564–8573. |
[22] | Barry M, Pearce H, Cross L, et al. (2016). Advances in nanotechnology for the treatment of osteoporosis. Curr Osteoporosis Rep 14: 87–94. |
[23] | Saini D, Fazil M, Ali MM, et al. (2014). Formulation, development and optimization of raloxifene-loaded chitosan nanoparticles for treatment of osteoporosis. Drug Deliv 7544: 1–14. |
[24] | Ponnapakkam T, Katikaneni R, Sakon J, et al. (2014) Treating osteoporosis by targeting parathyroid hormone to bone. Drug Discov Today 19: 204–208. |
[25] | Narayanan D, Anitha A, Jayakumar R, et al. (2013) In vitro and in vivo evaluation of osteoporosis therapeutic peptide PTH 1-34 loaded PEGylated chitosan nanoparticles. Mol Pharm 10: 4159–4167. |
[26] | Lindsay R, Krege JH, Marin F, et al. (2016) Teriparatide for osteoporosis: importance of the full course. Osteoporosis Int 2016: 1–16. |
[27] | Ong KL, Villarraga ML, Lau E, et al. (2010) Off-label use of bone morphogenetic proteins in the United States using administrative data. Spine J 35: 1794–1800. |
[28] | Carragee EJ, Hurwitz EL, Weiner BK (2011) A critical review of recombinant human bone morphogenetic protein-2 trials in spinal surgery: emerging safety concerns and lessons learned. Spine J 11: 471–491. |
[29] | Canalis E (2010) Update in new anabolic therapies for osteoporosis. J Clin Endocr Metab 95: 1496–1504. |
[30] | Tokatlian T, Segura T (2010) siRNA applications in nanomedicine. Wires Nanomed Nanobi 2: 305–315. |
[31] | Guo B, Wu H, Tang T, et al. (2012) A delivery system targeting bone formation surfaces to facilitate RNAi-based anabolic therapy. Nat Med 18: 307–314. |
[32] | Lin JH (1996) Bisphosphonates: a review of their pharmacokinetic properties. Bone 18: 75–85. |
[33] | Aoki K, Alles N, Soysa N, et al. (2012) Peptide-based delivery to bone. Adv Drug Deliver Rev 64: 1220–1238. |
[34] | Cenni E, Granchi D, Avnet S, et al. (2008) Biocompatibility of poly(d,l-lactide-co-glycolide) nanoparticles conjugated with alendronate. Biomaterials 29: 1400–1411. |
[35] | Choi SW, Kim JH (2007) Design of surface-modified poly(d,l-lactide-co-glycolide) nanoparticles for targeted drug delivery to bone. J Control Release 122: 24–30. |
[36] | Xinluan W, Yuxiao L, Helena NH, et al. (2015) Systemic drug delivery systems for bone tissue regeneration-a mini review. Curr Pharm Design 21: 1575–1583. |
[37] | Yokogawa K, Miya K, Sekido T, et al. (2001) Selective delivery of estradiol to bone by aspartic acid oligopeptide and its effects on ovariectomized mice. Endocrinology 142: 1228–1233. |
[38] | Pignatello R (2011) PLGA-alendronate conjugate as a new biomaterial to produce osteotropic drug nanocarriers. Biomaterials Applications for Nanomedicine. InTech. |
[39] | Heller DA, Levi Y, Pelet JM, et al. (2013) Modular "click-in-emulsion" bone-targeted nanogels. Adv Mater 25: 1449–1454. |
[40] | Willson TM, Henke BR, Momtahen TM, et al. (1996) Bone targeted drugs 2 synthesis of estrogens with hydroxyapatite affinity. Bioorg Med Chem Lett 6: 1047–1050. |
[41] | Yarbrough DK, Hagerman E, Eckert R, et al. (2010) Specific binding and mineralization of calcified surfaces by small peptides. Calcified Tissue Int 86: 58–66. |
[42] | Carrow JK, Gaharwar AK (2015) Bioinspired polymeric nanocomposites for regenerative medicine. Macromol Chem Phys 216: 248–264. |
[43] | Lee MS, Su CM, Yeh JC, et al. (2016) Synthesis of composite magnetic nanoparticles Fe3O4 with alendronate for osteoporosis treatment. Int J Nanomed 11: 4583–4594. |
[44] | Kedmi R, Ben-Arie N, Peer D (2010) The systemic toxicity of positively charged lipid nanoparticles and the role of Toll-like receptor 4 in immune activation. Biomaterials 31: 6867–6875. |
[45] | Hwang SJ, Lee JS, Ryu TK, et al. (2016) Alendronate-modified hydroxyapatite nanoparticles for bone-specific dual delivery of drug and bone mineral. Macromol Res 24: 623–628. |
[46] | Lu T, Ma Y, Hu H, et al. (2011) Ethinylestradiol liposome preparation and its effects on ovariectomized rats' osteoporosis. Drug Deliv 18: 468–477. |
[47] | Allen TM, Cullis PR (2013) Liposomal drug delivery systems: from concept to clinical applications. Adv Drug Deliver Rev 65: 36–48. |
[48] | Kundu SK, Sharma AR, Lee SS, et al. (2014) Recent trends of polymer mediated liposomal gene delivery system. Biomed Res Int 2014: 934605–934619. |
[49] | Madni A, Sarfraz M, Rehman M, et al. (2014) Liposomal drug delivery: a versatile platform for challenging clinical applications. Journal Pharm Pharm Sci 17: 401–426. |
[50] | Fang C, Shi B, Pei YY, et al. (2006) In vivo tumor targeting of tumor necrosis factor-α-loaded stealth nanoparticles: effect of MePEG molecular weight and particle size. Eur J Pharm Sci 27: 27–36. |
[51] | Zhang J, Li X, Huang L (2014) Non-viral nanocarriers for siRNA delivery in breast cancer. J Control Release 190: 440–450. |
[52] | Hughes J, Yadava P, Mesaros R (2010) Liposomal siRNA delivery. Meth Mol Biol 605: 445. |
[53] | Ropert C (1999) Liposomes as a gene delivery system. Braz J Med Biol Res 32: 163–169. |
[54] | Hengst V, Oussoren C, Kissel T, et al. (2007) Bone targeting potential of bisphosphonate-targeted liposomes. Preparation, characterization and hydroxyapatite binding in vitro. Int J Pharm 331: 224–227. |
[55] | Vert M, Schwach G, Engel R, et al. (1998) Something new in the field of PLA/GA bioresorbable polymers? J Control Release 53: 85–92. |
[56] | Mooney DT, Mazzoni CL, Breuer C, et al. (1996) Stabilized polyglycolic acid fibre based tubes for tissue engineering. Biomaterials 17: 115–124. |
[57] | Jiang T, Yu X, Carbone EJ, et al. (2014) Poly aspartic acid peptide-linked PLGA based nanoscale particles: potential for bone-targeting drug delivery applications. Int J Pharm 475: 547–557. |
[58] | Fu YC, Fu TF, Wang HJ, et al. (2014) Aspartic acid-based modified PLGA-PEG nanoparticles for bone targeting: in vitro and in vivo evaluation. Acta Biomater 10: 4583–4596. |
[59] | Cong Y, Quan C, Liu M, et al. (2015) Alendronate-decorated biodegradable polymeric micelles for potential bone-targeted delivery of vancomycin. J Biomat Sci Polym E 26: 629–643. |
[60] | Yi H, Wu LQ, Bentley WE, et al. (2005) Biofabrication with chitosan. Biomacromolecules 6: 2881–2894. |
[61] | Aktaş Y, Andrieux K, Alonso MJ, et al. (2005) Preparation and in vitro evaluation of chitosan nanoparticles containing a caspase inhibitor. Int J Pharm 298: 378–383. |
[62] | Schwarz K, Milne D (1972) Growth-promoting effects of Silicon in rats. Nature 239: 333–334. |
[63] | Jugdaohsingh R (2007) Silicon and bone health. J Nutr Health Aging 11: 99–110. |
[64] | Boyce BF, Yao Z, Xing L (2010) Functions of nuclear factor κB in bone. Ann Ny Acad Sci 1192: 367–375. |
[65] | Nanes MS (2003) Tumor necrosis factor-α: molecular and cellular mechanisms in skeletal pathology. Gene 321: 1–15. |
[66] | Yamaguchi M, Neale WM (2012) The intact strontium ranelate complex stimulates osteoblastogenesis and suppresses osteoclastogenesis by antagonizing NF-κB activation. Mol Cell Biochem 359: 399–407. |
[67] | Ha SW, Neale WM, Beck GR (2014) Bioactive silica nanoparticles promote osteoblast differentiation through stimulation of autophagy and direct association with LC3 and p62. ACS Nano 8: 5898–5910. |
[68] | Liu F, Fang F, Yuan H, et al. (2013) Suppression of autophagy by FIP200 deletion leads to osteopenia in mice through the inhibition of osteoblast terminal differentiation. J Bone Miner Res 28: 2414–2430. |
[69] | Baeza A, Manzano M, Colilla M, et al. (2016) Recent advances in mesoporous silica nanoparticles for antitumor therapy: our contribution. Biomater Sci 4: 803–813. |
[70] | Vallet RM, Balas F, Arcos D (2007) Mesoporous materials for drug delivery. Angew Chem 46: 7548–7558. |
[71] | Vallet RM, Rámila A, Del Real RP, et al. (2001) A new property of MCM-41: drug delivery system. Chem Mater 13: 308–311. |
[72] | Sun W, Han Y, Li Z, et al. (2016) Bone-targeted mesoporous silica nanocarrier anchored by zoledronate for cancer bone metastasis. Langmuir 32: 9237–9244. |
[73] | Baeza A, Colilla M, Vallet RM (2015) Advances in mesoporous silica nanoparticles for targeted stimuli-responsive drug delivery. Expert Opin Drug Del 12: 319–37. |
[74] | Paris JL, Torre PDL, Manzano M, et al. (2016) Decidua-derived mesenchymal stem cells as carriers of mesoporous silica nanoparticles. in vitro and in vivo evaluation on mammary tumors. Acta Biomater 33: 275–282. |
[75] | Jordan A, Scholz R, Maier HK, et al. (2006) The effect of thermotherapy using magnetic nanoparticles on rat malignant glioma. J Neuro-Oncol 78: 7–14. |
[76] | Figuerola A, Di Corato R, Manna L, et al. (2010) From iron oxide nanoparticles towards advanced iron-based inorganic materials designed for biomedical applications. Pharmacol Res 62: 126–143. |
[77] | Choi SY, Song MS, Ryu PD, et al. (2015) Gold nanoparticles promote osteogenic differentiation in human adipose-derived mesenchymal stem cells through the Wnt/beta-catenin signaling pathway. Int J Nanomed 10: 4383–4392. |
[78] | Lee NK, Choi YG, Baik JY, et al. (2005) A crucial role for reactive oxygen species in RANKL-induced osteoclast differentiation. Blood 106: 852–859. |
[79] | Sul OJ, Kim JC, Kyung TW, et al. (2010) Gold nanoparticles inhibited the receptor activator of nuclear factor-κb ligand (RANKL)-induced osteoclast formation by acting as an antioxidant. Biosci Biotech Bioch 74: 2209–2213. |
[80] | Lean JM, Jagger CJ, Kirstein B, et al. (2015) Hydrogen peroxide is essential for estrogen-deficiency bone loss and osteoclast formation. Endocrinology 146: 728–735. |
[81] | Sahana H, Khajuria DK, Razdan R, et al. (2013) Improvement in bone properties by using risedronate adsorbed hydroxyapatite novel nanoparticle based formulation in a rat model of osteoporosis. J Biomed Nanotechnol 9: 193–201. |
[82] | Lin L, Chow KL, Leng Y (2009) Study of hydroxyapatite osteoinductivity with an osteogenic differentiation of mesenchymal stem cells. J Biomed Mater Res Part A 89: 326–335. |
[83] | Webster TJ, Ergun C, Doremus RH, et al. (2000) Enhanced functions of osteoblasts on nanophase ceramics. Biomaterials 21: 1803–1810. |
[84] | Yamaguchi M (2010) Role of nutritional zinc in the prevention of osteoporosis. Mol Cell Biochem 338: 241–254. |
[85] | Ito A, Otsuka M, Kawamura H, et al. (2005) Zinc-containing tricalcium phosphate and related materials for promoting bone formation. Curr Appl Phys 5: 402–406. |
1. | Hai-yang Xu, Heng-you Lan, Fan Zhang, General semi-implicit approximations with errors for common fixed points of nonexpansive-type operators and applications to Stampacchia variational inequality, 2022, 41, 2238-3603, 10.1007/s40314-022-01890-7 | |
2. | Yan Wang, Rui Wu, Shanshan Gao, The Existence Theorems of Fractional Differential Equation and Fractional Differential Inclusion with Affine Periodic Boundary Value Conditions, 2023, 15, 2073-8994, 526, 10.3390/sym15020526 | |
3. | Sumati Kumari Panda, Thabet Abdeljawad, Fahd Jarad, Chaotic attractors and fixed point methods in piecewise fractional derivatives and multi-term fractional delay differential equations, 2023, 46, 22113797, 106313, 10.1016/j.rinp.2023.106313 | |
4. | Philip N. A. Akuka, Baba Seidu, C. S. Bornaa, Marko Gosak, Mathematical Analysis of COVID-19 Transmission Dynamics Model in Ghana with Double-Dose Vaccination and Quarantine, 2022, 2022, 1748-6718, 1, 10.1155/2022/7493087 | |
5. | SUMATI KUMARI PANDA, ABDON ATANGANA, THABET ABDELJAWAD, EXISTENCE RESULTS AND NUMERICAL STUDY ON NOVEL CORONAVIRUS 2019-NCOV/ SARS-COV-2 MODEL USING DIFFERENTIAL OPERATORS BASED ON THE GENERALIZED MITTAG-LEFFLER KERNEL AND FIXED POINTS, 2022, 30, 0218-348X, 10.1142/S0218348X22402149 | |
6. | Tahair Rasham, Separate families of fuzzy dominated nonlinear operators with applications, 2024, 70, 1598-5865, 4271, 10.1007/s12190-024-02133-0 | |
7. | Sumati Kumari Panda, Thabet Abdeljawad, A. M. Nagy, On uniform stability and numerical simulations of complex valued neural networks involving generalized Caputo fractional order, 2024, 14, 2045-2322, 10.1038/s41598-024-53670-4 | |
8. | Sumati Kumari Panda, Kumara Swamy Kalla, A.M. Nagy, Limmaka Priyanka, Numerical simulations and complex valued fractional order neural networks via (ɛ – μ)-uniformly contractive mappings, 2023, 173, 09600779, 113738, 10.1016/j.chaos.2023.113738 | |
9. | Hammed Anuoluwapo Abass, Maggie Aphane, Morufu Oyedunsi Olayiwola, An inertial method for solving systems of generalized mixed equilibrium and fixed point problems in reflexive Banach spaces, 2023, 16, 1793-5571, 10.1142/S1793557123502078 | |
10. | Yassine Adjabi, Fahd Jarad, Mokhtar Bouloudene, Sumati Kumari Panda, Revisiting generalized Caputo derivatives in the context of two-point boundary value problems with the p-Laplacian operator at resonance, 2023, 2023, 1687-2770, 10.1186/s13661-023-01751-0 | |
11. | Nihar Kumar Mahato, Bodigiri Sai Gopinadh, , 2024, Chapter 15, 978-981-99-9545-5, 339, 10.1007/978-981-99-9546-2_15 | |
12. | Sumati Kumari Panda, Velusamy Vijayakumar, Bodigiri Sai Gopinadh, Fahd Jarad, 2024, Chapter 6, 978-981-99-9545-5, 177, 10.1007/978-981-99-9546-2_6 | |
13. | Sumati Kumari Panda, Ilyas Khan, Vijayakumar Velusamy, Shafiullah Niazai, Enhancing automic and optimal control systems through graphical structures, 2024, 14, 2045-2322, 10.1038/s41598-024-53244-4 | |
14. | P K Santra, G S Mahapatra, Sanjoy Basu, Stability analysis of fractional epidemic model for two infected classes incorporating hospitalization impact, 2024, 99, 0031-8949, 065237, 10.1088/1402-4896/ad4692 | |
15. | Anupam Das, Bipan Hazarika, Mohsen Rabbani, 2024, Chapter 10, 978-981-99-9206-5, 165, 10.1007/978-981-99-9207-2_10 | |
16. | Tahair Rasham, Sumati Kumari Panda, Ghada Ali Basendwah, Aftab Hussain, Multivalued nonlinear dominated mappings on a closed ball and associated numerical illustrations with applications to nonlinear integral and fractional operators, 2024, 10, 24058440, e34078, 10.1016/j.heliyon.2024.e34078 | |
17. | Sumati Kumari Panda, Velusamy Vijayakumar, Kottakkaran Sooppy Nisar, Applying periodic and anti-periodic boundary conditions in existence results of fractional differential equations via nonlinear contractive mappings, 2023, 2023, 1687-2770, 10.1186/s13661-023-01778-3 | |
18. | Sumati Kumari Panda, Kottakkaran Sooppy Nisar, Velusamy Vijayakumar, Bipan Hazarika, Solving existence results in multi-term fractional differential equations via fixed points, 2023, 51, 22113797, 106612, 10.1016/j.rinp.2023.106612 | |
19. | Sumati Kumari Panda, Velusamy Vijayakumar, Results on finite time stability of various fractional order systems, 2023, 174, 09600779, 113906, 10.1016/j.chaos.2023.113906 | |
20. | Sumati Kumari Panda, A.M. Nagy, Velusamy Vijayakumar, Bipan Hazarika, Stability analysis for complex-valued neural networks with fractional order, 2023, 175, 09600779, 114045, 10.1016/j.chaos.2023.114045 | |
21. | Mokhtar Bouloudene, Fahd Jarad, Yassine Adjabi, Sumati Kumari Panda, Quasilinear Coupled System in the Frame of Nonsingular ABC-Derivatives with p-Laplacian Operator at Resonance, 2024, 23, 1575-5460, 10.1007/s12346-023-00902-z | |
22. | Sumati Kumari Panda, Vijayakumar Velusamy, Ilyas Khan, Shafiullah Niazai, Computation and convergence of fixed-point with an RLC-electric circuit model in an extended b-suprametric space, 2024, 14, 2045-2322, 10.1038/s41598-024-59859-x |
h | 0.03125 | 0.0625 | 0.125 | 0.25 | 0.5 | 1.0 |
Error for S(t) | 0.27673 | 0.55049 | 1.0889 | 2.1307 | 4.0817 | 7.5129 |
Error for I(t) | 0.14805 | 0.2971 | 0.59806 | 1.2115 | 2.4843 | 5.2145 |
Error for R(t) | 0.41436 | 0.83717 | 1.6766 | 3.3317 | 6.5556 | 12.717 |