A proposed fractional dynamic system and Monte Carlo-based back analysis for simulating the spreading profile of COVID-19.

This paper presents a dynamic system for estimating the spreading profile of COVID-19 in Thailand, taking into account the effects of vaccination and social distancing. For this purpose, a compartmental network is built in which the population is divided into nine mutually exclusive nodes, including susceptible, insusceptible, exposed, infected, vaccinated, recovered, quarantined, hospitalized, and dead. The weight of edges denotes the interaction between the nodes, modeled by a series of conversion rates. Next, the compartmental network and corresponding rates are incorporated into a system of fractional partial differential equations to define the model governing the problem concerned. The fractional degree corresponding to each compartment is considered the node weight in the proposed network. Next, a Monte Carlo-based optimization method is proposed to fit the fractional compartmental network to the actual COVID-19 data of Thailand collected from the World Health Organization. Further, a sensitivity analysis is conducted on the node weights, i.e., fractional orders, to reveal their effect on the accuracy of the fit and model predictions. The results show that the flexibility of the model to adapt to the observed data is markedly improved by lowering the order of the differential equations from unity to a fractional order. The final results show that, assuming the current pandemic situation, the number of infected, recovered, and dead cases in Thailand will, respectively, reach 4300, 4.5 × 10 6 , and 36,000 by the end of 2021.

Introduction

The mathematical modeling of infectious diseases has been the foundation of infectious disease epidemiology for more than a century [1]. In recent years, detailed scrutiny of such diseases has become extensive because of improvements in data analysis, computing methods, diagnostic tests, and genome sequencing. Advances in such methods enable researchers to deploy mathematical models to design practical schemes for disease control, and also develop and test mathematical and statistical methods [2-6]. Mathematical models that describe the spread of diseases are continually being developed and are playing a fundamental role in promoting public health strategies in many countries [7-10].

In December 2019, a novel coronavirus, officially called coronavirus disease 2019 (COVID-19) by the World Health Organization (WHO), led to an outbreak of atypical pneumonia, first in Wuhan, the capital city of Hubei Province in China, and then expeditiously flared out in the entire world [11]. As of October 01, 2021, there have been more than 230 million confirmed cases, together with nearly 4.8 million deaths reports in the world [12]. Throughout anti-epidemic struggles, apart from medical, biological, and scientific research, theoretical studies based on mathematical and statistical modeling may thus significantly contribute to the conception of outbreak attributes [13-15]. Forecasting the trajectory of the spread accordingly helps countries make informed decisions on the required actions to mitigate the spread of the disease.

Different models have thus far been proposed to analyze COVID-19, but compartmental methods that divide the population into compartments have dominated epidemiological modeling area [16-19]. These models assist in interpreting the way COVID-19 spreads and in forecasting regional pinnacles of the pandemic. Compartmental models are broadly based on studying the systems of differential equations. In this respect, differential equations concentrate on the rate of changes in a variable or a group of variables as time proceeds. The most common model is the SIR, which divides a population into susceptible (S), infected (I), and recovered (R). The popularity of this model stems from its simplicity in predicting a small number of parameters. There are a number of studies that predict the spread of COVID-19 by this method, e.g., see [20, 21], and [22].

Various models have been similarly extracted from the basic SIR model, which have thus far supplemented additional compartments to the SIR model [23]. One of the valuable models is the SEIR, namely, the susceptible-exposed-infected-recovered model, which has added the exposed compartment to the SIR structure [24]. The exposed group refers to the phase between the susceptible and the infected individuals. It comprises the population exposed to an infection but not yet infected. The SEIR model is thus capable of analyzing the spread of diseases more accurately than the basic SIR model [25]. Of note, a host of studies have analyzed the COVID-19 pandemic through modified versions of the SEIR model, e.g., see [26, 27] and [28].

Some studies have proposed the fractional SIR model, where one or, feasibly, both sub-compartments are raised by exponents that are usually less than unity [29]. The reason is that the infection transmits outwards from an infected population to the whole population through the early stages of an epidemic. Within this framework, where susceptible is much greater than infected it is better to scale them as a fractional power. Given the decreasing number of infected people over the general population in real world, the exponent is anticipated to be larger than 1/2. Furthermore, the power of the susceptible population, which is notably higher than that of the infected population, may be negligible at least in the initial stages of an epidemic.

Deploying fractional exponents stemmed from a growth model, called the Norton–Simon–Massagué (NSM) model [30]. This model was created to explain the growth of biological organisms through applying determined energy principles. The governing differential equation readswhere and measure anabolism growth and defuse, respectively. Equation 1 may be simplified as declaring the net growth rate of a biological organism outcome of the balance between synthetic and degradative mechanisms. Meanwhile, the rate is proportional to the growth volume G(t) with a power function and the rate of the latter process relates linearly with G(t). It is thus vital to mention that the two exceptional cases of Eq. 1 are power-law, where , and second-type growth, where , to explain the tumor growth [31]. Different types of fractional compartmental models have been applied to a number of infectious diseases, e.g., see [32-38].

Despite the fact that these models are convincing instruments and typically allow practical intuition into grasping the growth and the spread of diseases within designated populations, they suffer from one major flaw; that is, with the increase of the number of compartments, the number of model parameters increases, making it difficult for the model to properly fit the data. In other words, there is a trade-off between increasing the computational node and the regression accuracy. To cover new aspects of the problem, such as quarantine, hospitalization, and vaccination, it is necessary to add new components to the problem to make model more realistic [39, 40]. This, on the other hand, increases the number of model parameters, greatly adds to the complexity of the problem, and in turn reduces the accuracy of the fit.

To address this issue, this paper puts forward a time-variant SEIVR-QH network, comprised of nine nodes: (1) susceptible, (2) insusceptible, (3) vaccinated, (4) exposed, (5) infected, (6) quarantined, (7) recovered, (8) hospitalized and (9) dead. In a novel development, the differential equation is expressed in fractional form, described shortly. The conventional form of compartmental networks employs first-order differential equations. In contrast, introducing the fractional form here in this paper facilitates the formulation of differential equations in non-integer, fractional orders, which leads to significant flexibility for the model to describe the observed data. This model has several unknown parameters that must be estimated to fit the model to the actual data. It requires a multi-objective optimization to select the optimal parameters by minimizing the error between the actual data and the model outputs. In addition, this optimization process requires solving the system of fractional differential equations in each iteration which is a complex task. To this end, a two-step back analysis is proposed in this paper based on Monte Carlo sampling. A comparison is made here between the capability of the proposed method and that of other optimization techniques, such as the particle swarm optimization (PSO) used in [41]. First, PSO is a deterministic optimization algorithm and Monte Carlo sampling is an inherently probabilistic method. In other words, the system parameters in the latter are modeled as random variables. The back analysis proposed in this study is a two-step algorithm. That is, first, the Monte Carlo sampling method is used to identify the samples that lead to the closest realizations to the target time histories. The results are then introduced into the second stage of the back-analysis algorithm to obtain the optimal parameters. This leads to the probabilistic characteristics of the optimal model parameters, such as the mean and standard deviation, using the output of the first step. In addition, the deterministic optimal parameters can be obtained as the mean values. Therefore, the proposed back analysis can provide both probabilistic and deterministic parameters, while the PSO method only yields optimal deterministic results. Second, the modeling approach in the present paper is based on the fractional form of the system of differential equations. The performance of PSO in the optimization of fractional equations has not yet been evaluated. In contrast, Monte Carlo sampling method is robust in handling such complex systems of fractional networks proposed in this study. Finally, the fractional order of the proposed differential equations also needs to be optimized along with other model parameters. Hence, the problem here involves differential equations of varying order, which is intractable in the classical PSO. That is, the classical PSO does not allow the optimization of the order of differential equations as the change in the order will fundamentally change the objective function and affect the process of calculating the cost function. In contrast, Monte Carlo sampling is based on the generation of random samples and the fractional order is yet another random variable in this analysis, which can be readily handled.

The content of the paper is organized as follows. In Sect. 2, the proposed time-variant fractional SEIVR-QH network is presented, and the deterministic form of the problem is formulated. In Sect. 3, the main framework of the Monte Carlo-based back analysis is introduced and applied to the actual data of Thailand. In Sect. 4, a sensitivity analysis [42] is conducted to investigate the influence of the fractional order of the differential equation on the regression accuracy. Finally, the research is summarized and concluded in Sect. 5.

Deterministic fractional SEIVR-QH network

Fractional calculus

In this paper, the modeling formulation is derived by the Caputo derivative in fractional calculus. Several researchers have applied fractional derivatives in science and engineering, e.g., see [43].

Definition 1

Given a positive real number and an integrable function , the fractional integral of f of order is defined as [44]:where denotes the Gamma function.

Definition 2

The Caputo fractional derivative, where , is defined as [44]:

If f is a function of class , its fractional derivative is represented by the following expression:where represents the Gamma function.

Deterministic SEIVR-QH network

In this paper, an extended fractional SEIR network, called the SEIVR-QH, is applied to develop a deterministic model for the spread of COVID-19. In the definition of this network, the population is divided into nine nodes. First, the whole population is considered susceptible (S), which can become exposed (E) or can be immunized by observing social distancing (P) or being vaccinated (V) against the disease. Depending on the social exposure, exposed cases can be transmitted to the infected stage (I). If an individual is diagnosed with the disease, they will be quarantined (Q) based on the disease severity and infection spread. If the infection is severe, the patient will then be hospitalized (H). Consequent to treatment via hospital-based medical care or self-medication, nearly all quarantined cases recover (R) and some die (D) because of the disease severity. The relationship between these nine states is depicted in the network shown in Fig. 1.

Fig. 1

Relationship between different nodes

The mathematical relationships among the nine states in the model are expressed by a system of fractional differential equations, as follows:where is the protection rate, indicates the infection rate, stands for the vaccination rate, is the average latency time, denotes the rate at which infected cases enter into quarantine, refers to hospitalization rate, is the time-variant quarantine recovery rate, is the time-variant mortality quarantine rate, and and are hospitalization recovery and mortality rates, respectively. The constraint , , is applied, where N is the entire target population.

Theorem 1

There is a unique solution for the fractional SEIVR-QH model. Moreover, the solution is non-negative.

Proof

Let and for , f(x) is the Lebesgue measure. Then,with , . If , then the function f is non-decreasing, and if , then the function f is non-increasing for all . Afterwards, all the conditions in Theorem 3.1 from [45] are met and there is a unique and non-negative function in , which is the solution of Eq. (1). There is also a need to show that the domain is positively invariant. Since the following holds:on each hyper-plane bounding the non-negative orthant, the vector field points into .

Theorem 2

Consider the initial value problem ofEqs. (5)–(14), with satisfying , to evaluate its equilibrium points. Let

Then, the equilibrium points are

A disease-free equilibrium

An endemic equilibrium point , where , , and .

Numerical solution

The system of fractional differential equations is written in the matrix form, as follows:To solve the matrix form of the given system, the Adams–Bashforth–Moulton predictor–corrector method is used in two steps. First, the prediction step calculates a rough approximation of the desired quantity, typically employing an explicit method. Second, the corrector step refines the initial approximation recruiting another means, typically an implicit method. For more details, see [46].

Monte Carlo sampling for back analysis

At this point, the formulation of SEIVR-QH network, the fractional form of the differential equations, and its solution are established. The problem at hand is to estimate the unknown weights of the network, namely , , , , , , , , , , , and , such that the proposed SEIVR-QH network best describes the observed data. To this end, a novel two-step optimization algorithm is established in this section using Monte Carlo-based back analysis to obtain the optimal weights when is assumed as unity [47-50]. A sensitivity analysis is then used in the next section to investigate the influence of fractional order in the regression accuracy. To implement the back-analysis approach, the parameters involved in the problem, i.e., weights of the network, are defined as random variables. This makes it possible to consider a range of variation for these parameters and evaluate the predictions obtained by their various combinations. To better understand the parameters’ combination, suppose the parameters are placed in a lattice, in which the x-axis and y-axis are the parameters within the arranged set of ranges, then the variety of these parameters with different values is supposed as a random walk on this lattice. The range of variations defined for each random variable is shown in Table 1.

Table 1

Random variables and their corresponding ranges of variation

No.	Parameter	Distribution function	Minimum	Maximum
1	β	Uniform	0.0	10.0
2	α	Uniform	0.0	0.5
3	γ	Uniform	0.0	0.05
4	δ	Uniform	0.0	0.5
5	λ_0	Uniform	0.0	0.5
6	λ_1	Uniform	0.0	0.5
7	κ_0	Uniform	0.0	0.5
8	κ_1	Uniform	0.0	0.5
9	ζ	Uniform	0.0	0.1
10	ρ	Uniform	0.0	0.1
11	ψ	Uniform	0.0	0.1
12	ϕ	Uniform	0.0	0.1

Next, 10,000 realizations are generated for each random variable. This process takes 30 min on a PC equipped with a 64 GB of random access memory (RAM) and an Intel Core i7-5820K central processing unit (CPU). Then, by introducing each set of realizations into the system of differential equations, a new SEIVR-QH problem is established that is solved using the method describe previously. This leads to 10,000 distinct predictions of the spread of the disease over time. The resulting database of I, R, D, and V is then used as the input of the back-analysis algorithm in the next step. The process described for implementing the Monte Carlo sampling is schematically shown in Fig. 2.

Fig. 2

Steps of using Monte Carlo sampling in combination with the SEIVR-QH network

To implement the back-analysis method, the actual observations of COVID-19, i.e., I, R, and D in Thailand, are compiled from the World Health Organization (WHO) online database [12]. In addition, the number of daily vaccinations is obtained from the Johns Hopkins University database [51]. The whole dataset is available on “Mendeley Data” [52]. The results are plotted in Fig. 3 from 2021-02-28 to 2021-08-30. This figure also demonstrates intervals of above and below the observed data. This is here considered as an allowable interval, and used as a filter to select the closest cases to the observed data among the 10,000 realizations of I, R, D, and V. This means that the allowable area is defined by a sample-selection criterion that leads the algorithm towards finding the best available fit. The process described for the selection/rejection of samples is schematically shown in Fig. 4.

Fig. 3

The allowable interval defined above and below the observed data for a infected, b recovered, c dead, and d vaccinated cases

Fig. 4

Schematic diagram of the proposed back-analysis method

By applying the proposed filter to the data, 40 out of 10,000 samples fall within the allowable interval. Next, a new criterion is defined for a second-round selection to find the best sample from the previously selected 40 samples. For this purpose, the root mean square error (RMSE) between the predicted time series of I, R, D, and V in each sample and the observed time series is calculated for each of the 40 samples, as follows:where , , , and are the RMSE of the sample for the infected, recovered, dead, and vaccinated, respectively; I(t), R(t), D(t) and V(t) are the actual measurements of infected, recovered, dead, and vaccinated cases, respectively; and , , , and are the time series of infected, recovered, dead, and vaccinated cases for the sample, respectively, computed using the proposed SEIVR-QH model. The results for , , , and are shown in columns 2–5 and 8–11 of Table 2. Next, by combining the RMSEs of I, R, D, and V using a square root of the sum of squares, the final selection criterion, denoted by , is defined:The results of are shown in columns 6 and 12 of Table 2. The previously selected 40 samples are then sorted according to and the sample corresponding to the lowest , i.e., sample 4 in Table 2, is selected as the best fit. The optimal values of the parameters correspond to this sample are reported in Table 3. The time series of I, R, D, and V corresponding to the optimal parameters are also shown in Fig. 5.

Table 2

RMSE for 40 samples

No.	θ_1 (×10^4)	θ_2 (×10^4)	θ_3 (×10^6)	θ_4 (×10^2)	θ_t (×10^6)	No.	θ_1 (×10^4)	θ_2 (×10^4)	θ_3 (×10^6)	θ_4 (×10^2)	θ_t (×10^6)
1	3.96	2.05	4.69	3.09	4.69	21	9.45	1.87	4.35	2.86	4.35
2	5.07	2.26	4.63	2.90	4.63	22	3.09	2.38	4.60	3.49	4.60
3	10.30	2.31	4.39	2.77	4.40	23	8.36	2.47	3.85	3.25	3.85
4	2.50	1.83	3.52	1.39	3.52	24	5.60	2.44	4.16	3.50	4.16
5	5.15	2.02	4.25	3.20	4.25	25	10.70	2.26	4.44	2.74	4.45
6	5.19	1.99	4.65	2.84	4.65	26	8.33	2.08	4.59	2.36	4.59
7	6.77	2.45	4.55	2.03	4.55	27	2.51	2.38	4.07	3.28	4.07
8	4.11	2.42	4.58	2.70	4.58	28	3.86	2.51	4.17	2.28	4.17
9	4.79	2.00	4.44	3.22	4.44	29	7.62	2.52	4.09	3.10	4.09
10	7.83	2.44	3.83	2.47	3.83	30	9.81	2.29	4.54	2.53	4.54
11	4.27	2.45	4.23	2.94	4.23	31	6.65	2.31	4.67	1.95	4.67
12	10.00	2.15	4.44	2.64	4.44	32	7.08	2.53	4.17	1.82	4.17
13	8.13	2.02	4.66	2.36	4.66	33	7.84	2.52	4.18	1.95	4.18
14	2.53	2.34	4.45	2.79	4.45	34	3.48	2.46	4.22	2.59	4.22
15	4.17	2.39	4.38	2.42	4.38	35	6.65	2.35	3.90	3.44	3.90
16	9.64	2.07	4.52	3.19	4.52	36	4.50	2.25	4.56	3.06	4.56
17	5.36	2.28	4.61	3.08	4.61	37	9.45	2.33	4.48	3.47	4.48
18	5.41	2.11	4.46	1.69	4.46	38	9.44	2.06	4.73	2.63	4.73
19	7.52	2.38	4.61	3.34	4.61	39	10.00	2.42	4.62	3.03	4.62
20	4.36	2.31	4.17	3.19	4.17	40	5.81	2.13	4.54	1.93	4.54

Table 3

The optimal values of parameters

No.	Parameter	Value	No.	Parameter	Value
1	β	5.1769	7	κ_0	0.0360
2	α	0.1599	8	κ_1	0.1239
3	γ	0.0294	9	ζ	0.0061
4	δ	0.3246	10	ρ	0.0837
5	λ_0	0.1773	11	ψ	0.0010
6	λ_1	0.3034	12	ϕ	0.0015

Fig. 5

Model predictions corresponding to the optimal parameters

Sensitivity analysis

As described in Sect. 2, the fractional form of the SEIVR-QH model, i.e., , provides a new degree of freedom to the compartmental system and provides the algorithm with further flexibility compared to the conventional differential equations of the first order, i.e., . This subsequently enables the system to better fit the actual data. Since the back analysis presented in this paper is implemented by four objective functions that must be optimized simultaneously to ensure that the number of infected, recovered, dead, and vaccinated cases matches the actual data, this flexibility is key in achieving accurate results. To make use of the fractional differential equation system to fit the SEIVR-QH network to the collected data of COVID-19 in this paper, the order of the equations is defined as a decision variable that assumes values less than unity. Thereafter, the new fractional problem is solved and the resulting fit is assessed for improvement in the prediction accuracy. To make use of the fractional differential equation system to fit the SEIVR-QH network to the collected data of COVID-19, the order of the equations is defined as a decision variable that assumes values less than unity. Then, the new fractional problem is solved and the resulting fit is assessed for improvement in the prediction accuracy.

For this purpose, assuming the optimal parameters obtained in Sect. 3, 11 different values of the decision variable that represent the fractional order are considered, including 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95 and 1. Next, a fractional form of SEIVR-QH problem is defined corresponding to each and analyzed separately by the method presented in Sect. 2. The results of each case are then stored in the form of time series I, R, and D. Table 4 shows the , , and , as well as for different values of . In addition, the corresponding time series of I, R, and D are plotted against different in Fig. 6.

Table 4

RMSE values for different

No.	η	θ_1 (×10^4)	θ_2 (×10^4)	θ_3 (×10^6)	θ_4 (×10^2)	θ_t (×10^6)
1	0.50	4.20	2.16	4.07	1.90	4.07
2	0.55	5.29	2.29	3.89	2.19	3.89
3	0.60	5.55	2.12	3.54	1.77	3.54
4	0.65	2.96	2.12	3.59	1.41	3.59
5	0.70	3.14	2.26	3.87	2.19	3.87
6	0.75	5.72	2.21	4.10	2.00	4.10
7	0.80	4.49	2.22	4.20	2.22	4.20
8	0.85	3.39	2.27	3.43	2.03	3.43
9	0.90	2.12	1.81	2.58	1.15	2.58
10	0.95	4.74	2.21	3.37	1.90	3.37
11	1.00	2.50	1.83	3.52	1.39	3.52

Fig. 6

The resulting fit on the data for different

As seen, the best result is obtained for . In other words, assuming , a better fit to data is obtained and is minimized. Compared to the assumption of , the value of has changed from to , i.e., a reduction of 27% in the error measured. According to the new results, the number of infected, recovered, and dead at the beginning of 2022 in the Thailand will, respectively, reach about 4300, and 36,000. In addition, by November 2022, the number of infected will fall below 1500, indicating that the disease will be under control in the Thailand. It is noteworthy that these results are based on the fit obtained by the fractional model, which is more accurate and reliable than the results reported in Sect. 3.

Conclusion

In this paper, a compartmental network is established using fractional SEIVR-QH to predict the spread profile of COVID-19 in Thailand. To this end, the proposed network is employed in a Monte Carlo-based back analysis in which network weights, including , , , , , , , , , , , and , are defined as random variables. By generating a set of random samples for each random variable, the Monte Carlo-back analysis is used to find realizations of parameters whose predictions fall within a desired interval around the observed data. Thereafter, the RMSE criterion is applied to the selected samples to find the best fit to the Thailand data. Finally, a sensitivity analysis is conducted to investigate the effect of the order of the fractional operator on the accuracy of the fit. The primary results of this study are as follows:Although the compartmental networks benefit from simple configuration and straightforward implementation, they suffer from a common problem; that is, if compartments and computational nodes are added, their complexity is greatly increased, making it difficult to optimize the network to adapt to the data. The method proposed in this paper, however, enables the compartmental network to be expanded to the required size without concerning the complexity of the network and the size of the differential equation system. In other words, one can simply add the required nodes, such as reinfection after recovery and vaccine side effects, and then use Monte Carlo sampling method to optimize the network. This will be studied by authors in the future. Moreover, SEIVR-QH parameters are assumed as variables in the back analysis and incorporated in matching of the time histories. However, in addition to the model parameters, the fractional order of the differential equations is also defined as random variables and employed in the formulation of the Monte Carlo sampling method. This allows for a higher flexibility to fit the model to the data and thus, increases the model accuracy. Secondly, the probabilistic characteristics of the order of the differential equation, namely the mean and the standard deviation, are obtained from the optimal results of the first stage of the back analysis. This research is a practical step in defining the randomized fractional problem in solving the differential equations of disease spread. This is the subject of ongoing research by the authors on modeling the spread of COVID-19.

Monte Carlo-back analysis identified 40 samples in close proximity of the observed data, having RMSE values between and . The best fit between them provides the optimal network weights as infection rate, , protection rate, , infectious to quarantine rate, , average latent time, , time-dependent recovery rate, and , time-dependent mortality rate, and , hospitalization rate, , vaccination rate, , hospitalization to recovery rate, , and hospitalization to mortality rate, .

The results of sensitivity analysis showed that the accuracy of data fitting is greatly improved by changing the fractional order, so that the of the best fractional fit, i.e., , reaches about , while that for the ordinary case, i.e., , is bound to . This means that the fractional form of the SEIVR-QH network improves the accuracy by about 27%.

The optimal output of the fractional SEIVR-QH network showed that by the end of 2021, the number of people infected (active), recovered and dead in Thailand will reach 4300, and 36,000, respectively.

The resulting number of infected cases also shows that the number of active cases fall below 1500 by November 2022, indicating that the disease will be under control in Thailand.