Citation: Masaki Sekiguchi, Emiko Ishiwata, Yukihiko Nakata. Dynamics of an ultra-discrete SIR epidemic model with time delay[J]. Mathematical Biosciences and Engineering, 2018, 15(3): 653-666. doi: 10.3934/mbe.2018029
[1] | Yoichi Enatsu, Yukihiko Nakata, Yoshiaki Muroya . Global stability for a class of discrete SIR epidemic models. Mathematical Biosciences and Engineering, 2010, 7(2): 347-361. doi: 10.3934/mbe.2010.7.347 |
[2] | Ke Guo, Wanbiao Ma . Global dynamics of an SI epidemic model with nonlinear incidence rate, feedback controls and time delays. Mathematical Biosciences and Engineering, 2021, 18(1): 643-672. doi: 10.3934/mbe.2021035 |
[3] | Yukihiko Nakata, Ryosuke Omori . The change of susceptibility following infection can induce failure to predict outbreak potential by $\mathcal{R}_{0}$. Mathematical Biosciences and Engineering, 2019, 16(2): 813-830. doi: 10.3934/mbe.2019038 |
[4] | A. Q. Khan, M. Tasneem, M. B. Almatrafi . Discrete-time COVID-19 epidemic model with bifurcation and control. Mathematical Biosciences and Engineering, 2022, 19(2): 1944-1969. doi: 10.3934/mbe.2022092 |
[5] | Edward J. Allen . Derivation and computation of discrete-delayand continuous-delay SDEs in mathematical biology. Mathematical Biosciences and Engineering, 2014, 11(3): 403-425. doi: 10.3934/mbe.2014.11.403 |
[6] | Zhen Jin, Zhien Ma . The stability of an SIR epidemic model with time delays. Mathematical Biosciences and Engineering, 2006, 3(1): 101-109. doi: 10.3934/mbe.2006.3.101 |
[7] | Jianquan Li, Zhien Ma, Fred Brauer . Global analysis of discrete-time SI and SIS epidemic models. Mathematical Biosciences and Engineering, 2007, 4(4): 699-710. doi: 10.3934/mbe.2007.4.699 |
[8] | Zenab Alrikaby, Xia Liu, Tonghua Zhang, Federico Frascoli . Stability and Hopf bifurcation analysis for a Lac operon model with nonlinear degradation rate and time delay. Mathematical Biosciences and Engineering, 2019, 16(4): 1729-1749. doi: 10.3934/mbe.2019083 |
[9] | Oren Barnea, Rami Yaari, Guy Katriel, Lewi Stone . Modelling seasonal influenza in Israel. Mathematical Biosciences and Engineering, 2011, 8(2): 561-573. doi: 10.3934/mbe.2011.8.561 |
[10] | Shiqiang Feng, Dapeng Gao . Existence of traveling wave solutions for a delayed nonlocal dispersal SIR epidemic model with the critical wave speed. Mathematical Biosciences and Engineering, 2021, 18(6): 9357-9380. doi: 10.3934/mbe.2021460 |
Differential equations are used for modeling phenomena in many disciplines e.g. economy, physics, engineering, biology and chemistry. Since exact (analytical) solutions are not available in many cases if the mathematical model is given as a nonlinear system, numerical solutions have enhanced our understanding of the mathematical model [7].Traditional numerical schemes such as Euler's method and Runge-Kutta method, however, may induce numerical instability, as the discretization would change properties of solutions, such as stability and positivity, of the original model, thus numerical scheme should be carefully chosen to preserve nature of the original system, see e.g. [4,8].
The authors in [10] propose a system of difference equations as a discrete counterpart of a continuous epidemic model, describing disease transmission dynamics in continuous time. Qualitative properties of the model such as positivity, boundedness and convergence of the solutions are investigated. As a continuation work of [10] in [11] a general system of difference equations is analyzed. In both papers, proving convergence of the solution, when a parameter called the basic reproduction number is greater than one, seems to be a challenging problem, while it is known that the corresponding continuous model has an equilibrium that is globally stable. A "good" discretization is found in [23,24,5] for a class of epidemic models formulated by delay differential equations in [2,3,29]. The authors in [24] prove uniform persistence of the solution, corresponding to a result in [29] for a continuous SIRS model. The authors in [5] prove global stability of the endemic equilibrium by a Lyapunov function, corresponding to the result established in [14] (see also [6]). The discretization used in [24,5] is a variation of backward Euler's method and is indeed known as Micken's nonstandard finite difference method [15]. See also [12,16] for the application of Micken's nonstandard finite difference method to ordinary differential equation models. It is also known that some discrete-time epidemic models exhibit periodic and chaotic behavior [1].
Ultra-discretization is proposed as a procedure to obtain a discrete system, where unknown variables also take discrete values, thus a cellular automaton is defined, see [21,19]. In [27,22,20] the authors study discrete and ultra-discrete models for epidemic models. In those papers it is shown that simple ultra-discrete models can capture the disease transmission dynamics, which is seen in the original differential equation models. In [27] the authors find conserved quantities for some special cases. Cellular automata have been used to model disease transmission dynamics, see e.g. [25,26] and references therein. Since cellular automata are computational models, in general, it is not straightforward to perform a mathematical analysis, in order to provide theoretical insights into simulation studies. Our analytical approach in this paper for the ultra-discrete model may be used to complement numerical simulation studies for some cellular automaton models.
In this paper we start with a special case of the model considered in [5,24] for a discrete analogue of an epidemic model formulated by a system of delay differential equations in [2,3,29,14]. As delay differential equation form an infinite dimensional dynamical system [9], after discretizing the system, we obtain a system of difference equations of higher order, which is slightly complicated compared to the model considered in [27]. From such a model we derive an ultra-discrete model, which is formulated as a couple of piecewise linear difference equations (see e.g. [13]). Due to time delay, integrability can not be usually expected and the application of the ultra-discretization to a non-integrable system is still limited, but see also [17,18]. We here prove that the ultra-discrete model exactly has the same threshold property regarding global attractivity as in [5,14]. In Section 4 we visualize the convergence of the solution in a two-dimensional lattice and observe some interesting convergence patterns. For a special initial condition, we derive a simple recurrence relation for the solution in Lemma 3. The relation can explain the illustrated convergence pattern.
The paper is organized as follows. In Section 2 we introduce a system of difference equations, which is a special case of the model studied in [5,24] for discrete analogue of continuous models. Applying a variable transformation together with taking a limit, we derive an ultra-discrete model. In Section 3 global behavior of the solutions is discussed. We prove that the model exhibits the threshold behavior, similar to the difference equation studied in [5]. We here find that a subsequence of the solution has a monotone property and this monotonicity is used for the proof. In Section 4 we illustrate the solution behavior in a two-dimensional lattice. To explain the convergence pattern, we consider a special initial condition and derive a simple recurrence relation. We summarize our results in Section 5.
Our starting point is the following system of difference equations
Sn+1−Sn=M−MSn+1−BSn+1In−ω, | (1a) |
In+1−In=BSn+1In−ω−MIn+1−ΓIn+1, | (1b) |
Rn+1−Rn=ΓIn+1−MRn+1, | (1c) |
with a positive initial condition. System (1) is a special case of the model considered in [24,5]. See also [6] for a model with nonlinear incidence. Here
R0=BM+Γ. |
One can prove that system (1) always has the disease free equilibrium given by
(M+ΓB,MB(BM+Γ−1),ΓB(BM+Γ−1)) |
if
Theorem 2.1.
(ⅰ) If
(ⅱ) If
System (1) is proposed as a discrete analogue of the following disease transmission dynamics model in continuous time:
S′(t)=M−MS(t)−BS(t)I(t−τ), | (2a) |
I′(t)=BS(t)I(t−τ)−MI(t)−ΓI(t), | (2b) |
R′(t)=ΓI(t)−MR(t), | (2c) |
where
Let us now derive an ultra-discrete model from system (1), following the same procedure as in [27,22,20]. Since the third equation of (1) does not appear in the first and second equations for
Sn=exn/ϵ and In=eyn/ϵ, |
and parameters
M=eμ/ϵ,B=eβ/ϵ and Γ=eγ/ϵ. |
Notice that (1a) and (1b) are equivalently written as
Sn+1=Sn+M1+M+BIn−ω, | (3a) |
In+1=In+BSn+1In−ω1+M+Γ. | (3b) |
Applying the variable transformation to (3) with letting
xn+1=max(μ,xn)−max(0,μ,β+yn−ω), | (4a) |
yn+1=max(yn,β+xn+1+yn−ω)−max(0,μ,γ). | (4b) |
The key relation used here is the following limit
limϵ→+0ϵlog(eA/ϵ+eB/ϵ)=max(A,B) |
for
The initial condition of the system (4) is given as
x0 and y−j∈R for j∈{0,1,⋯,ω}. | (5) |
Moreover, parameters
xn+1=max(μ,xn)−max(μ,β+yn−ω), | (6a) |
yn+1=max(yn,β+xn+1+yn−ω)−γ. | (6b) |
In this section, we study the global asymptotic behavior of the solutions of (6).
Lemma 3.1. For any solutions, there exists
Proof.
Let us assume that
xn+1=xn−max(μ,β+yn−ω)≤xn−μ, |
thus
xk+1=μ−max(μ,β+yk−ω)≤0. |
Thus we get that
xn+1=μ−max(μ,β+yn−ω), |
from (6a)
β+xn+1+yn−ω=μ+β+yn−ω−max(μ,β+yn−ω)≤μ |
follows. Therefore, from (6b), we have
yn+1≤max(yn−γ,μ−γ). | (7) |
If
From Lemma 3.1, without loss of generality, we can set the initial condition as
x0≤0 and y−j≤μ−γ≤0, j∈{0,1,⋯,ω}. |
Note that Lemma 3.1 implies
{(x,y)∈R2 | x≤0,y≤μ−γ≤0} | (8) |
is an invariant set. To discuss global attractivity of equilibria, we now consider the following system
xn+1=μ−max(μ,β+yn−ω), | (9a) |
yn+1=max(yn,β+xn+1+yn−ω)−γ | (9b) |
in (8).
Theorem 3.2. Let us assume that
xn=0forn≥1andlimn→∞yn=−∞. |
Proof.
Since for any
max(μ,β+yn−ω)=μ. |
Therefore it follows that
yn+1=max(yn−γ,β−γ+yn−ω). | (10) |
Let
Ym:=max0≤j≤ωym(ω+1)−j for m∈N+. |
Note that
y(m+1)(ω+1)−j≤Ym+max(−γ,β−γ). |
From (10), we have
y(m+1)(ω+1)−ω=ym(ω+1)+1=max(ym(ω+1)−γ,β−γ+ym(ω+1)−ω)≤max(Ym−γ,β−γ+Ym)=Ym+max(−γ,β−γ). |
For some
y(m+1)(ω+1)−j≤Ym+max(−γ,β−γ). | (11) |
Then using (10) and (11) we obtain
y(m+1)(ω+1)−j+1=max(y(m+1)(ω+1)−j−γ,β−γ+ym(ω+1)−j+1)≤Ym+max(−γ,β−γ). |
By mathematical induction, it holds that
Ym+1≤Ym+max(−γ,β−γ). |
Now it is obvious that
If
Proposition 1. Let us assume that
Proof.
Let
μ−max(μ,β+μ−γ)=μ−(β+μ−γ)=−β+γ,max(μ−γ,β−β+γ+μ−γ)−γ=max(μ−γ,μ)−γ=μ−γ. |
Proposition 2. Let us assume that
yn+1={μ−γif yn−ω≥μ−βmax(yn−γ,β−γ+yn−ω)if yn−ω<μ−β. |
Proof.
Assume that
yn+1=max(yn−γ,μ−γ)=μ−γ. |
On the other hand, assume that
We now show that every solution converges to the non-trivial equilibrium.
Theorem 3.3 Let us assume that
limn→∞xn=−β+γandlimn→∞yn=μ−γ. |
Proof . Let
yℓ:=(yℓ(ω+1),yℓ(ω+1)−1,⋯,yℓ(ω+1)−ω) |
for
yℓ(ω+1)−k≥β−γ+y(ℓ−1)(ω+1)−k |
for
limℓ→∞yℓ(ω+1)−k=μ−γ,k∈{0,1,2,⋯,ω}, |
i.e. each component of
xn+1=μ−(β+μ−γ)=−β+γ. |
Corresponding to the first and second parts of Theorem 2.1, we show the threshold behavior in Theorems 3.2 and 3.3 for the ultra-discrete epidemic model (9).
In Figures 1 and 2 we plot
In Figure 2, we set
We here visualize the convergence of a solution using a two-dimensional lattice. We consider the case that
u−jm:=ym(ω+1)−j for m∈N+ and j∈{0,1,⋯,ω} |
to set a sequence with two variables
u−jm={μ−γif u−jm−1≥μ−βmax(u−j−1m−γ,u−jm−1+β−γ)if u−jm−1<μ−β, | (12) |
for
u−ωm={μ−γif u−ωm−1≥μ−βmax(u0m−1−γ,u−ωm−1+β−γ)if u−ωm−1<μ−β. | (13) |
In Figure 3, specifying parameters and the initial condition, we write a numerical value of
To explain the pattern we consider a special solution and derive an explicit recurrence relation for the solution in Proposition 3. For
u−j0=p for j∈{0,1,⋯,ω}∖{k},u−k0=q, | (14) |
with
p<q<μ−β. | (15) |
Lemma 4.1. Assume that
p+kγ<q. | (16) |
Then it holds
u−j1=p+β−γ for j∈{k+1,k+2,⋯,ω} | (17a) |
u−j1=q+β−γ−(k−j)γ for j∈{0,1,⋯,k}. | (17b) |
Proof.
Firstly, we show (17a) by mathematical induction. From (14), (15) and (16), one has
u−ω1=max(p−γ,p+β−γ)=p+β−γ. |
Suppose that
u−j1=p+β−γ for some j∈{k+2,k+3,⋯,ω}. |
From (12) one can see
u−j+11=max(p+β−2γ,p+β−γ)=p+β−γ. |
Therefore we get (17a).
Next we show (17b) by mathematical induction. From (14) one can see
u−k1=max(u−k−11−γ,u−k0+β−γ)=max(p+β−γ−γ,q+β−γ)=q+β−γ. |
Next suppose that
u−j1=q+β−γ−(k−j)γ for some j∈{1,2,⋯,k}. |
Then from (12), one can see
u−j+11=max(q+β−γ−(k−j+1)γ,p+β−γ)=q+(β−γ)−(k−j+1)γ. |
We note that
q+(β−γ)−(k−j+1)γ−(p+β−γ)=q−p−(k−j+1)γ≥q−p−kγ>0 |
holds by (16). Therefore, we get (17b).
Proposition 3. Assume that (16) and
q<p+β−γ+(k+1)γ | (18) |
hold. Then it holds that
u−jm=β−γ+u−jm−1 if u−jm−1<μ−β, | (19) |
for
Proof.
Note that (19) with
u−j2=β−γ+u−j1. | (20) |
Here from (14) one has that
u−ω2=max(u01−γ,u−ω1+β−γ)=u−ω1+β−γ. |
Here using (17a), (17b) and (18), we know that
u−ω1+β−γ−(u01−γ)=p−q+β−γ+(k+1)γ>0. | (21) |
Next suppose that
u−j2=u−j1+β−γ for some j∈{1,2,⋯,ω}. | (22) |
From (12) and (22), one has
u−j+12=max(u−j2−γ,u−j+11+β−γ)=max(u−j1+β−γ−γ,u−j+11+β−γ)=max(u−j1−γ,u−j+11)+(β−γ). |
Let us show that
u−j+11−(u−j1−γ)≥0 for j∈{1,2,⋯,ω}, | (23) |
so that
u−j+12=β−γ+u−j+11 | (24) |
holds. For
u−j1−γ=p+(β−γ)−γ,u−j+11=q+(β−γ). |
Thus we have
Finally, assuming
u−jm=u−jm−1+β−γ for j∈{0,1,⋯,ω} | (25) |
and
u−ωm+1=max(u0m−γ,u−ωm+β−γ). |
Note that from (25) and (21), one can see
u−ωm+β−γ−(u0m−γ)=u−ω1+β−γ−(u01−γ)>0. |
Therefore we get
u−ωm+1=u−ωm+β−γ. |
For some
u−jm+1=u−jm+β−γ | (26) |
holds. From (12) one has
u−j+1m+1=max(u−jm+1−γ,u−j+1m+β−γ). |
Note that from (26) and (25) one can see
u−jm+1−γ=u−j1+m(β−γ)−γ,u−j+1m+β−γ=u−j+11+m(β−γ). |
Then from (23) one can see
u−j+1m+β−γ−(u−jm+1−γ)=u−j+11−(u−j1−γ)≥0. |
Therefore we get (19) for
When
u−j+11=u−j1−γ<u−j1 for j∈{1,2,⋯,k}. |
In this way
Although here we consider a special solution such that the initial condition is given as in (14), we can show that linear combination of two solutions can produce a complicated pattern for the convergence of the solution. Consider two solutions
w−jm:=max(u−jm,v−jm) | (27) |
is shown to be a solution of (12) and (13) as follows. Assume that
wjm=max{max(uj−1m−γ,ujm−1+β−γ),max(vj−1m−γ,vjm−1+β−γ)}=max(uj−1m−γ,ujm−1+β−γ,vj−1m−γ,vjm−1+β−γ)=max(max(uj−1m,vj−1m)−γ,max(ujm−1,vjm−1)+β−γ)=max(wj−1m−γ,wjm−1+β−γ). |
Assume that either
Let the solution depicted in Figure 3(a) be
In this paper we consider an ultra-discrete model with time delay. The model is derived from a discrete epidemic model studied in [5,24]. In Theorems 3.2 and 3.3, we show that the ultra-discrete model also has the threshold property concerning global attractivity of equilibria, similar to the discrete epidemic model [5] and the continuous epidemic model [14]. For the proof of global attractivity of the non-trivial equilibrium in Theorem 3.3, we reduce the system (9) to the scalar difference equation in Proposition 2 and then use a certain monotone property of the solution, which is our important finding.
In Section 4 we further derive a simple recurrence relation for the solution, assuming a special condition for the initial condition. The relation derived in Proposition 3 clearly shows that each component monotonically increases towards the equilibrium. Two-dimensional lattice seems to be an informative tool to illustrate such a convergence patten.
For the ultra-discrete model (9), a linear combination of two solutions is shown to be a solution. It would be interesting to investigate the structure of the solutions of ultra-discrete models. Figure 3 may remind of an elementary cellular automaton (see e.g. rule 252 in [28]). Exploring a possible connection to elementary cellular automaton is our future work.
We are very grateful to the anonymous referee for the evaluation of our paper and for the constructive critics.
The second author was supported by JSPS Grant-in-Aid for Scientific Research (C) JP26400212. The third author was supported by JSPS Grant-in-Aid for Young Scientists (B) 16K20976.
[1] | [ L. J. S. Allen, Some discrete-time SI, SIR and SIS epidemic models, Math. Bio., 124 (1994): 83-105. |
[2] | [ E. Beretta,Y. Takeuchi, Global stability of an SIR epidemic model with time delays, J. Math. Biol., 33 (1995): 250-260. |
[3] | [ E. Beretta,T. Hara,W. Ma,Y. Takeuchi, Global asymptotic stability of an SIR epidemic model with distributed time delay, Nonl. Anal., 47 (2001): 4107-4115. |
[4] | [ R. M. Corless,C. Essex,M. A. H. Nerenberg, Numerical methods can suppress chaos, Phys. Lett. A, 157 (1991): 27-36. |
[5] | [ Y. Enatsu,Y. Nakata,Y. Muroya, Global stability for a class of discrete SIR epidemic models, Math. Bio. and Eng., 7 (2010): 347-361. |
[6] | [ Y. Enatsu,Y. Nakata,Y. Muroya,G Izzo,A Vecchio, Global dynamics of difference equations for SIR epidemic models with a class of nonlinear incidence rates, J. Diff. Equ. Appl., 18 (2012): 1163-1181. |
[7] | [ S. Elaydi, An Introduction to Difference Equations, Springer-Verlag, New York, 2005. |
[8] | [ D. F. Griffiths,P. K. Sweby,H. C. Yee, On spurious asymptotic numerical solutions of explicit Runge-Kutta methods, IMA J. Numer. Anal., 12 (1992): 319-338. |
[9] | [ J. K. Hale and S. M. Verduyn Lunel, Introduction to Functional Differential Equations, Applied Mathematical Sciences Vol 99, Springer, 1993. |
[10] | [ G. Izzo,A. Vecchio, A discrete time version for models of population dynamics in the presence of an infection, J. Comput. Appl. Math., 210 (2007): 210-221. |
[11] | [ G. Izzo, Y. Muroya and A. Vecchio, A general discrete time model of population dynamics in the presence of an infection Disc. Dyn. Nat. Soc. , (2009), Art. ID 143019, 15pp. |
[12] | [ L. Jódar,R. J. Villanueva,A. J. Arenas,G. C. González, Nonstandard numerical methods for a mathematical model for influenza disease, Math. Comput. Simul., 79 (2008): 622-633. |
[13] | [ C. M. Kent, Piecewise-defined difference equations: Open Problem, 'Bridging Mathematics, Statistics, Engineering and Technology, 55-71, Springer Proc. Math. Stat., 24, Springer, New York, 2012. |
[14] | [ C. C. McCluskey, Complete global stability for an SIR epidemic model with delay -distributed or discrete time delays, Nonl. Anal. RWA, 11 (2010): 55-59. |
[15] | [ R. E. Mickens, Discretizations of nonlinear differential equations using explicit nonstandard methods, J. Comput. Appl. Math., 110 (1999): 181-185. |
[16] | [ S. M. Moghadas,M. E. Alexander,B. D. Corbett,A. B. Gumel, A positivity-preserving Mickens-type discretization of an epidemic model, J. Diff. Equ. Appl., 9 (2003): 1037-1051. |
[17] | [ K. Matsuya,M. Murata, Spatial pattern of discrete and ultradiscrete Gray-Scott model, DCDS-B, 20 (2015): 173-187. |
[18] | [ K. Matsuya and M. Kanai, Exact solution of a delay difference equation modeling traffic flow and their ultra-discrete limit, arXiv: 1509.07861 [nlin. CG]. |
[19] | [ K. Nishinari,D. Takahashi, Analytical properties of ultradiscrete Burgers equation and rule-184 cellular automaton, J. Phys. A, 31 (1998): 5439-5450. |
[20] | [ A. Ramani,A. S. Carstea,R. Willox,B. Grammaticos, Oscillating epidemics: A discrete-time model, Phys. A, 333 (2004): 278-292. |
[21] | [ T. Tokihiro,D. Takahashi,J. Matsukidaira,J. Satsuma, From soliton equations to integrable cellular automata through a limiting procedure, Phys. Rev. Lett., 76 (1996): 3247-3250. |
[22] | [ J. Satsuma,R. Willox,A. Ramani,B. Grammaticos,A. S. Carstea, Extending the SIR epidemic model, Phys. A, 336 (2004): 369-375. |
[23] | [ M. Sekiguchi, Permanence of some discrete epidemic models, Int. J. Biomath., 2 (2009): 443-461. |
[24] | [ M. Sekiguchi,E. Ishiwata, Global dynamics of a discretized SIRS epidemic model with time delay, J. Math. Anal. Appl., 371 (2010): 195-202. |
[25] | [ G. C. Sirakoulis,I. Karafyllidis,A. Thanailakis, A cellular automaton model for the effects of population movement and vaccination on epidemic propagation, Ecol. Mod., 133 (2000): 209-223. |
[26] | [ S. H. White,A. Martin del Rey,G. Rodríguez Sánchez, Modeling epidemics using cellular automata, Appl. Math. Comp., 186 (2007): 193-202. |
[27] | [ R. Willox,B. Grammaticos,A. S. Carstea,A. Ramani, Epidemic dynamics: Discrete-time and cellular automaton models, Phys. A, 328 (2003): 13-22. |
[28] | [ S. Wolfram, Statistical mechanics of cellular automata, Rev. Mod. Phys., 55 (1983): 601-644. |
[29] | [ T. Zhang,Z. Teng, Global behavior and permanence of SIRS epidemic model with time delay, Nonl. Anal. RWA., 9 (2008): 1409-1424. |
1. | Fernando Córdova-Lepe, Rodrigo Gutiérrez, Karina Vilches-Ponce, Analysis of two discrete forms of the classic continuous SIR epidemiological model, 2020, 26, 1023-6198, 1, 10.1080/10236198.2019.1696323 | |
2. | P. S. Knopov, A. S. Korkhin, Dynamic models of epidemiology in discrete time taking into account processes with lag, 2023, 2195-268X, 10.1007/s40435-023-01135-3 | |
3. | Yue Zhang, Xue Li, Xianghua Zhang, Guisheng Yin, Lei Chen, Stability and Hopf Bifurcation Analysis of an Epidemic Model with Time Delay, 2021, 2021, 1748-6718, 1, 10.1155/2021/1895764 |