A comprehensive analysis of COVID-19 nonlinear mathematical model by incorporating the environment and social distancing

This research conducts a detailed analysis of a nonlinear mathematical model representing COVID-19, incorporating both environmental factors and social distancing measures. It thoroughly analyzes the model’s equilibrium points, computes the basic reproductive rate, and evaluates the stability of the model at disease-free and endemic equilibrium states, both locally and globally. Additionally, sensitivity analysis is carried out. The study develops a sophisticated stability theory, primarily focusing on the characteristics of the Volterra–Lyapunov (V–L) matrices method. To understand the dynamic behavior of COVID-19, numerical simulations are essential. For this purpose, the study employs a robust numerical technique known as the non-standard finite difference (NSFD) method, introduced by Mickens. Various results are visually presented through graphical representations across different parameter values to illustrate the impact of environmental factors and social distancing measures.


Model building process
Upon a comprehensive review of the research paper as cited above, we came to the conclusion that their suggested model was not adequately scrutinized and needs further investigation and analysis.We refined the model by reverting the recovered class into the susceptible class due to considerations regarding the weakened immune system, and subsequently made additional adjustments to the model outlined in system (1).We then revisit the updated model given in (2) and establish a detailed analysis for the stability and numerical results by using NSFD method. (1) + ϕE − µS, − (ϕ + µ + υ)E, Table 1.Summarization of the key results of each literature.
Year Key data sources highlighting our research findings 2021  Chien and Shateyi 30 , primarily employed global stability analysis using V-L matrices to investigate the transmission dynamics of Babesiosis 2021 Shao and Shateyi 49 , employed mathematical modeling, specifically utilizing Lyapunov functions and V-L matrices, to analyze the global stability of the SEIRS epidemic model with a non-linear incident rate, and validated the findings through numerical simulations 2020 Mwalili et al. 1 , utilized mathematical modeling, particularly a modified SEIR compartmental model, alongside numerical methods like fourth and fifth order Runge-Kutta techniques, to forecast COVID-19 epidemic patterns and analyzed the basic reproduction number using the next generation matrix approach 2020 Masoumnezhad et al. 50, studied global stability analysis by utilizing a modification of the V-L matrix method, which combines Lyapunov functions and V-L matrices to analyze the mathematical model of an infectious disease 2018 Verma and Kayenat 32 , analyzed and established convergence of Mickens' type non-standard finite difference (NSFD) schemes, applied them to Lane-Emden equations, and compared the numerical results with existing analytical solutions and standard finite difference (FD) schemes, showing NSFD's ability to efficiently approximate solutions near singular points where FD fails 2017 Zahedi and Kargar 48 , used a methodologies involved mathematical modeling of HIV/AIDS transmission dynamics, incorporating constant controls, V-L stable matrices, and numerical simulations to analyze global stability 2012 Mickens and Washington 36 , created an NSFD scheme for an SIRS model of respiratory virus transmission, preserving conservation laws and ensuring positive solutions for all time step sizes 2012 Liao and 28 , employed mathematical epidemiological models integrating human and environmental factors, using V-L stable matrices alongside classical Lyapunov functions to analyze global stability, yielding new results for three-and four-dimensional systems, validated via numerical simulations www.nature.com/scientificreports/

Compartmental structure
The population is compartmentalized based on disease status and the pathogen population into susceptible (S), exposed (E), symptomatic infectious ( I S ), asymptomatic infectious ( I A ), recovered (R), and pathogen (P) categories, defined in Table 2.A graphical representation illustrating the proposed model governing the system of differential equations outlined in (2) is shown in Fig. 1.Each compartment is governed by a differential equation describing the rate of change of individuals or the pathogen population over time.Susceptible individuals can transition to exposed through contact with infected individuals or the pathogen, influenced by transmission dynamics.Exposed individuals progress to infectious states, symptomatic or asymptomatic, at specific rates.Recovery and disease-induced mortality rates determine transitions from infectious to recovered states.Pathogen population dynamics are affected by transmission from infectious individuals and decay processes.
The model integrates parameters governing disease transmission, progression rates, recovery rates, mortality rates, and pathogen dynamics nomenclature is given in Table 2.This model construction process allows for a comprehensive representation of infectious disease dynamics, facilitating insights into transmission patterns, intervention strategies, and disease control measures.Table 2. Physical illustration of nomenclature of the model (2).

Model testing approach
To ensure the reliability and applicability of our proposed model, we conducted a rigorous testing process, employing various methodologies to assess its predictive capabilities and robustness.The testing framework comprised the following: The model's predictive accuracy has been validated by comparing its simulations with past epidemiological data obtained from 1 .Through ongoing adjustments to the model's parameters, we ensured its close alignment with observed trends in disease incidence, prevalence, and other relevant metrics.This process confirmed the reliability of our model in forecasting epidemiological outcomes.
A comprehensive sensitivity analysis has been conducted to evaluate the impact of parameter variations on model outcomes.By systematically varying key parameters within plausible ranges, we assessed the sensitivity of the model to different inputs and identified parameters with the most significant influence on model predictions.This analysis enabled us to better understand the model's behavior under varying conditions and to identify potential sources of uncertainty.
Finally, to ensure the robustness and credibility of our model.The feedback received from peer reviewers strengthened the validity and reliability of our model and its implications for public health policy and practice.

Feasibility of solution
After describing the human population in system (2), it is imperative to demonstrate that, for all t ≥ 0 , the state parameters S, E, I A , I S , R, andP are positive.For every t ≥ 0 , the solution with positive beginning data remains positive and is bounded.Now from system (2), clearly one has that dN dt = b − µN − ((δ − τ S ) + (δ − τ A ) + µ P P) and sup t−→+∞ N ≤ b µ .The feasible region is described as At this point, (3) is positively invariant with respect to (2).This indicates that all solutions of the system with (S, E, I A , I S , R, P) ∈ R 6 + remain in and thus the suggested model ( 2) is epidemiologically well formulated.DFE Solving equations given in system (2) we get the DFE points by setting S ′ , E ′ , I ′ A , I ′ S , R ′ and P ′ all equal to zero.In this case I ′ A = I ′ S = P ′ = 0 , which implies that E = 0 and R = 0 too.Hence, we get the DFE E 0 = ( b µ , 0, 0, 0, 0, 0).EE To deduce EE point, let use S ′ , E ′ , I ′ A , I ′ S , R ′ and P ′ in system (2) in right sides and consider From 3rd and 4th equations in the system (4), we get and (2) Adding first two equations of system (4) while using R * = ωE * , we get

Computation of R 0
The basic reproduction rate, often represented as R 0 , quantifies the average count of subsequent infections initi- ated by one person within a population that is entirely susceptible.This metric serves as a gauge to determine if the infection is likely to propagate within the population.The calculation of R 0 employs the methodology known as the "next generation matrix approach".Since the first and fifth equations in system (2) are independent of each other, we have consequently reduced system (2) to the following: Let x = (E, I A , I S , P) T , then the model can be written as d dt (x) = F (x) − V(x) , where www.nature.com/scientificreports/Suppose F and V represent the Jacobian matrices of F and V at DFE respectively, then we obtain FV −1 as where c 1 = µ + υ , c 2 = µ + δ + γ S + τ S , and c 3 = µ + δ + γ A + τ A .The reproduction number R 0 is the largest eigenvalue of FV −1 , given by The notation stands for basic reproductive number for human R H 0 and for pathogens R P 0 .Our remarks are recorded as: and therefore Here, R 0 is composed of two portions, which correspond to the two COVID-19 transmission modes.In addition, we present 3D profile of R 0 corresponding to various parameters to investigate the trends of parameters involve in (13) in Fig. 2 as follows:

Local stability of DFE and EE
Theorem 3.1 The DFE is locally asymptotically stable for R 0 < 1.

Proof
The jacobian matrix with respect to the system (2) is given by The characteristic polynomial of J 0 about (S 0 , E 0 , I 0 S , I 0 A , R 0 , P 0 ) is given by Here with It follows from ( 14) that 1 = −µ , 2 = −(µ + ϕ) and 3 = −µ P are the eigenvalues of J 0 , its other eigenval- ues are the roots of (15).Furthermore, it is clear that both a > 0 and b > 0 .Considering the high sensitivity of parameters γ S , τ S , and β 2 , it is important to note that (γ to maintain the local asymptotic stability of the DFE provided that R 0 < 1 .Failure to meet these conditions would render the DFE locally asymptotically unstable provided that R 0 > 1 .Consequently, we conclude that c > 0 .This analysis establishes that ab − c > 0 .This indicates that the real parts of the roots of (15) are negative.Thus, the system's single positive equilibrium point of the system (2) is locally asymptotically stable according to the Hurwitz criterion.This implies that the DFE is locally asymptotically stable.This completes the proof.

Local stability analysis of EE at E *
In this subsection, we establish the asymptotic stability of ( 2) at E * .

Theorem 3.2
The EE denoted by E * is locally asymptotically stable with R 0 > 1.
Proof At E * , the jacobian matrix of (2) can be deduced as www.nature.com/scientificreports/For simplicit y, let The matrix shown above can be expressed in the following manner: here in this case We see that one eigenvalue of Jacobian matrix J * is clearly negative, i.e. 1 = −µ P < 0 .The remaining eigenval- ues can be tested by observing that det(J * ) and Tr(J * ) have different signs.Now from the above jacobian matrix J * , we have Considering the high sensitivity of parameters η , c 2 , and c 3 , it is important to note that η > c 2 c 2 −c 3 (assuming that c 2 = c 3 ) to maintain the local asymptotic stability of the EE provided that R 0 > 1 .Failure to meet these condi- tions would render the EE locally asymptotically unstable.In this context, we see that det(J * ) and Tr(J * ) have opposite signs which insures that the eigenvalues of J * have negative real parts.This implies that the system (2) is locally asymptotically stable at E * .

Global stability of DFE and EE
In this section, we investigate the V-L stability of matrix A defined in (16) to determine the global stability of the endemic equilibrium E * .The following definitions and preliminary lemmas are necessary prerequisites for this examination.Definition 4.1 30,[48][49][50] .Consider a square matrix, denoted as M, endowed with the property of symmetry and positive (negative) definiteness.In this context, the matrix M is succinctly expressed as M > 0 or ( M < 0).Definition 4.2 30,[48][49][50] .We write a matrix A n×n > 0(A n×n < 0) if A n×n is symmetric positive (negative) definite.Definition 4.3 30,[48][49][50] .If a positive diagonal matrix H n×n exists, such that HA + A T H T < 0 then a non-singular matrix A n×n is V-L stable.Definition 4.4 30,[48][49][50] .If a positive diagonal matrix H n×n exists, such that HA + A T H T < 0 (> 0) then, a non singular matrix A n×n is diagonal stable.

Global stability of DFE
We employ the same methodology used in 51 to obtain the necessary outcomes for globally asymptotically.

Proof Taking the Lyapunov function corresponding to all dependent classes of the model
Taking derivative with respect to t yields After utilizing the values derived from system (2) and performing requisite calculations, we obtain: Using in the above equation gives Thus, the system (2) is globally asymptotically stable at E 0 , whenever R 0 < 1 .

Global stability analysis of the EE at E *
We propose the following Lyapunov function in order to prove the global stability at E * .
Here, m 1 , m 2 , m 3 , w 4 , m 5 , and w 6 are non-negative constants.By taking the derivative of L with respect to time along the trajectories of the model ( 8), one has Here It is important to note that matrix A (n−1)×(n−1) is derived from matrix A by removing its final row and column.
In accordance with the works of Zahedi and Kargar 48 , Shao and Shateyi 49 , Masoumnezhad et al. 50, Chien and Shateyi 30 , and Parsaei et al. 29 , we introduce a series of lemmas and theorems to examine the global stability of the EE at E * .
In order to make the calculation easier, we present the matrix ( 16) in more simplified form as given below: H e r e , a = It is important to note that all the parameters a, b, c, d, e, f, g, h, i, j, k, l, m, n, p, and q are positive.Throughout this paper, it is important to note that a − e and −a + e are always equal to µ and −µ respectively.Moreover, we conclude that all the diagonal elements (the negative values) of the matrix A are negative, which is a good sign for stability.We need to compute the eigenvalues to make a definitive conclusion for the global stability of E * .(      Lemma 4. 8 The matrix E given by C −1 in Lemma 4.4 exhibits diagonal stability.
Proof It is obvious that E 44 > 0 , by applying Lemma 4.2, our task is to establish the diagonal stability of the reduced matrix I = E and its inverse I = E −1 .This demonstration serves as the completion of the proof.Thus the matrix I and I can be obtained from the matrix E given in Lemma 4.4, we have Here and Since det(E) > 0 , and i 33 > 0 .Therefore I 33 > 0. (17), and if (µ , then the following conditions holds: . This implies that From ( 19) and ( 20), we get From ( 18) and ( 21), we deduced that (iv) From (i) and (iii), one can easily see that fhj > b(hi
Lemma 4. 15 The matrix I given by E in Lemma 4.9 exhibits diagonal stability.
Vol.:(0123456789)In summary of the preceding discussions, we've drawn these conclusions regarding the global asymptotic stability of the EE.

Sensitivity analysis
Here in this section, we conduct sensitivity analysis to explore how the model responds to variations in parameter values, aiming to understand which parameters significantly influence the transmission of the disease, specifically the reproductive number ( R 0 ).To perform this analysis, we employ the normalized forward sensitivity index method, as outlined in 52 .This method involves determining the ratio of the relative change in a variable to the relative change in the corresponding parameter.Alternatively, sensitivity indices can be derived using partial derivatives when the variable is a differentiable function of the parameter.
Therefore, for our proposed model, we present the normalized forward sensitivity index of " R 0 " using the given role:

Numerical results and discussion
We simulate our results by extending the NSFD method.The mentioned scheme is convergent, consisted and have many dynamical properties (see 32 ).Mickens's 34 established a general NSFD method for a coupled system.The reason behind of develing such schemes was that to avoid inaccuracies in the standard methods like Euler, RKM methods.Here, consider a problem 35 described by we define t i+1 = t 0 + nh , where h is positive step size, let the discretization of y at t n be y n ≈ y(t n ) , then the the problem (25) becomes The procedure ( 26) is called NSFD method, because • In the discretized version of ( 26), the usual denominator h is replaced by a nonnegative function say φ(h); such that • The nonlinear function θ is approximated by nonlocal way.
The major advantages of the aforesaid method are given as (see details 33 ): • The NSFD method avoids numerical instabilities even on using very small step size we need to implement the scheme.• First, the NSFD method behaves more or equally well with an order RK4 technique most of the time.
Although the numerical algorithmic complexity is similar to that of the Euler scheme, there is a significant computational benefit.• Second, a system that provides good dynamic behavior even for huge time increments is produced by main- taining the dynamical restrictions.Additionally, this offers a significant computational and implementation advantage.
To deduce the NSFD algorithm for our model (2).We can write in differences form S ′ as Using ( 27) in ( 2), we get such that ζ(h) = 1 − exp(−h) (see 36 ), and h → 0 , then ζ(h) → 0 .After rearranging the terms, we get the algo- rithm in simplified form as From the numerical scheme we get (30), the positivity condition of the NSFD method clearly can be observed.We use the said method and the following numerical values given in Table 3 for the simulations to present graphically our results. ,

Case II
When α 2 = 0.05, 0.10, 0.15 , we describe the population dynamics of various classes using that values of α 1 = 0.10 in Figs. 10, 11, 12, 13, 14, 15.As time increases in the populations from 0 to 100 days, Fig. 10 shows this.As illustrated in Fig. 11, the number of susceptible individuals rapidly drops during the first ten days, while the number of exposed individuals rapidly increases during the first twenty days as a result of interaction with infected individuals.For the first 20 days after infection, the number of infected individuals is also growing quite quickly.Similarly, both symptomatic and asymptomatic infected classes are rapidly growing (see Figs. 12 and 13).As seen in Fig. 15, there is also a rise in pathogen concentration.As seen in Fig. 14, the population of the recovered class is likewise growing steadily.Here, for the first 25 days, the spread of diseases increases as the infectious environment ratio does.On reducing the social interaction with infectious individual, the chances of getting infection is also rising.As at α 2 = 0.05 , there is high risk to get infection.But on increasing the values of α 2 = 0.10 , and α 2 = 0.15 , the social distance increases the chances of getting infection will decrease.

Conclusion
The present article conducts a comprehensive review aimed at understanding the global dynamics of COVID-19.
The researchers identified equilibrium points, calculated the basic reproductive rate, and assessed the model's stability at disease-free and endemic equilibrium states, both locally and globally.The main focus lies in examining COVID-19 dynamics using a suggested model, with particular attention to the characteristics of V-L matrices.A series of lemmas and theorems based on V-L matrices theory are presented in this regard to study the global stability of the endemic equilibrium.Non-standard finite difference techniques are employed for numerical analysis, contributing to ongoing efforts to comprehend the pandemic's intricacies.The paper also elaborates on the effects of social distance and interaction with infected environments.Precautionary measures such as increasing social distances, wearing face masks, sanitation, and avoiding large gatherings are highlighted as effective strategies for reducing the risk of infection.These findings contribute to the broader understanding of COVID-19 dynamics and provide insights into potential strategies for mitigating its spread.As future work, the authors plan to utilize the Volterra-Lyapunov (V-L) matrix theory approach to investigate the global stability of other infectious disease models that have not yet been thoroughly studied using the said methodology.

Theorem 4 . 6
The matrix A 6×6 defined in(17) is V-L stable.Proof Obviously −A 66 > 0 .Consider C = − A be a 5 × 5 matrix that is produced by eliminating the final row and column −A .Utilizing Lemma 4.2, we can demonstrate the diagonal stability of matrices C = − A and C = − A −1 .Let examine the matrix C = − A and C = − A −1 , from(17)

Theorem 4 . 7
When the value of R 0 is greater than 1, the endemic equilibrium E * = (S * , E * , I * S , I * A , R * , P * ) in model (2) exhibits global asymptotic stability within .Proof Building upon above Theorem and a sequence of lemmas, it is assured that the EE of model (2) achieves global asymptotic stability.

Figure 3 .
Figure 3.The population dynamics of the proposed model (2).

Figure 4 .Figure 5 .
Figure 4. Dynamics of the susceptible class of the proposed model (2).
The dynamic of the refined model governs in the form of differential equations given in(2).
But finding the eigenvalues of the matrix A is quite time consuming.Therefore, we introduce a series of lemmas and theorems given below to examine whether all the eigenvalues of the matrix A have negative real parts or not?.