In this paper, we use the method of stochastization of one-step processes to demonstrate the advantages of a stochastic representation of the system. To achieve this goal, the high-level scientific computing language Julia was chosen. Because the chosen approach is based on the paradigm of analytical-numerical calculations, then its implementation will be performed using auxiliary packages that are domain-specific extensions (DSL): Catalyst.jl and DiferentialEquations.jl. As a result of this work, the method of stochastization of one-step processes was applied to the SIR epidemic model. mathematic modeling, diferential equations, schemes of kinetics reactions, stochastization of one-step Information and Telecommunication Technologies and Mathematical Modeling of High-Tech Systems (ITTMM-2021), 0000-0003-1202-2787 (A. V. Fedorov); 0000-0002-3036-0117 (A. O. Masolova); 0000-0001-7141-7610 (A. V. Korolkova); 0000-0002-0877-7063 (D. S. Kulyabov) Workshop Proceedings htp:/ceur-ws.org CEUR Workshop Proceedings (CEUR-WS.org) ISN1613-073
The ordinary diferential equations describe the deterministic behavior of the observed model which seldom does not represent the actual behavior of the system. Based on the previous research we suppose that stochastic diferential equations can give better results.
One of the ways of obtaining SDE describing a system is to build self-consistent one-step processes. This method is based on the systems of interaction that allow for the opportunity to obtain SDE. This method is divided into two approaches: analytical (scheme of interaction, program code) and numerical (program code).
The goal of the research is the automatization of obtaining the equations and their solutions. For this purpose, the current task is solved in one computing environment.
Workshop on information technology and scientific computing in the framework of the XI International Conference https://yamadharma.github.io/ (D. S. Kulyabov) CEUR
The solution of this task is based on the paradigm of the analytical-numerical computations. For its realization, we chose the programming language Julia [1] which supports the creation of the subject-oriented plug-ins (DSL) [2]. Additionally, we are using the Catalyst.jl package which consists of the analytical and numerical modules, where the numerical one is based on the DiferentialEquations.jl plug-in.
The following algorithm had been developed: the entry point is the building of a set of kinetic reactions (based on the considered method of stochastization of one-step processes), describing the intercommunication of objects inside the model, then this set is interchanged into a scheme of interaction. Further, we conduct analytical computations, which interchange the obtained scheme into the SDE system. The SDE system is solved with the usage of the DiferentialEquations.jl packaging. As a result, we obtain graphs for the estimation of quality and quantity behavior of the system’s characteristics.
This complex helps solve the problem of applying the analytical-numerical systems for obtaining the solution of our problems. Although it lacks the ability to obtain the SDE with a self-consistent stochastic part. The Catalyst.jl package needs upgrading by including program code in it to generate a self-consistent stochastic pat of the SDE. 2. Stochastization of the SIR model In chemical kinetics, the reaction of interaction of substances can be represented as a timedependent process. Equations are one of the ways to represent the dependences of the concentrations of substances. In this case, the diferential equations can be represented in the form of kinetic reaction schemes. A complex of simple chemical reactions can be a complex reaction that is its composition.
There is a small amount of work in which the analytical approach is used in conjunction with the Julia language and the Catalyst.jl helper library used to model equations describing complex chemical reactions [3]. where () is the number of susceptible individuals to the disease at the time , () is the number of infected individuals at the time , () is the number of afected individuals at the time , = () + () + ()
, is the coeficient of the intensity of interaction between healthy individuals and sick individuals (which subsequently leads to a disease), is the coeficient of the intensity of recovery of sick individuals.
On the example of the SIR epidemic model (
= ,
,
− ,
(
Based on the first equation of the system, we can say that the decrease in the number of susceptible individuals is directly proportional to the number of interactions with sick individuals. Those. a healthy individual becomes infected after direct contact.
From the second equation it follows that a decrease in the number of sick individuals is directly proportional to an increase in the number of recovered individuals. It is also worth noting that there is a similar relationship regarding the rate of disease of individuals due to the number of their contacts with susceptible ones.
And finally, the third equation shows a direct relationship between sick and recovered (i.e., the number of recovered is equal to the number of infected),
Summing up the description of the model, we can say that the process of spreading the epidemic consists of two stages: 1. Healthy individuals that are susceptible turn to ill; 2. Sick individuals who are was infected turn to recover.
+ −→ 2 ,
−→ .
(
Now, for the model under consideration, we should write down the interaction scheme (
The first step is to introduce the vector which is responsible for describing the states of the
system, in our case = (, , )
. Next, we introduce the operators , (where i is the reaction
number), denoting the initial and final states of the system, respectively, as well as the operator
, obtained from their diference. Such operators are calculated using the following algorithm:
1. Compiled interaction diagram describing the flow of specific reactions in the system;
2. For the first reaction from the interaction scheme, the operators и
1
1 will look like
this (
1
(
To find these operators, write down the coeficient of the presence of parameters on the left side of the reaction. The form of writing the sought operators is based on the introduced vector of system states = (, , )
, where the presence of each of the indicated
parameters matters. The parameters and are present before the reaction occurs in a
single copy, so we get the operator 1 = (
The next step in the implementation of this technique is to calculate the intensities of the
transitions + , − (
interaction scheme, the underlined index is the tensor component аnd +/− — the direction of the reaction («+» — forward reaction, «−» — backward reaction):
The following are formulas for calculating the intensities of transitions (
! ∏ =1 ( − )! − , = +
! ∏ =1 ( − )! .
Find the transition rates to our model, they are calculated separately for each of the reaction
schemes of cooperation. Let us write down the intensity of the transition + 1 (
(
! ∗ ( − 1)! ( − 1)! = ∗ ∗ ( − 1)! ( − 1)! ∗ ∗ ( − 1)! ( − 1)! = .
The inverse intensity − 1 transition does not exist, since there is no backlash for our model.
Similarly, we write out + 2 (
Obtaining the master kinetic equation (MKE) is an optional step, so let’s move on to obtaining stochastic equations. Within the framework of this technique, they can be written in the Fokker-Planck form, which corresponds to the Langevin form. Obtaining the equation in the Fokker-Planck form is also an intermediate step; for brevity, in this work we will not write it out and write out the equation in the Langevin form right away. To do this, we need to calculate the drift vector and the difusion matrix .
General formulas (
= (
2
=1
=
∑ +
∗
=
+
1 ∗ 1 = ∗ (
Based on the obtained drift vector (11) and the difusion matrix (12) we obtain an equation in the Langevin form (14): (13) (14)
() = () + () , ( ) = ( − ) + (− + − ) ( 1 2) . 3
Thus, the process of stochastization of the equations of the SIR epidemic model was performed. 2.2. Obtaining solutions of the system of stochastic equations The next step in this work will be to obtain a solution for the system of stochastic equations and to visualize them: −1 0 0 1 0 1. To do this, we will set additional system parameters: beta = 0.1 gamma = 0.01 N = sum(u0) #disease rate #recovery rate Thus, a value for the disease rate was set equal to 0.1. For the recovery rate value value was set to 0.05. Coficient N is the sum of all the main parameters of the system — , , . 2. Write the system of equations in the program form: function sir(du, u, beta, gamma, N, t) du[1] = -((beta*u[1]*u[2])/N) du[2] = ((beta*u[1]*u[2])/N) - (gamma*u[2]) du[3] = gamma*u[2] end In the numerical approach, implemented using the DiferentialEquations.jl library, the model equations are specified as a function. It receives the following parameters as input: • du - array of technical parameters responsible for diferentiation; • u - array of basic parameters described in the section ??; • beta, gamma, N - coeficients described above; • t - time parameter; 3. Similarly, we declare a function that represents the stochastic part of the model: function sir_noise(du, u, beta, gamma, N, t) du[1] = sqrt(abs( ((beta*u[1]*u[2])/N)* (gamma*u[2]) )) du[2] = sqrt(abs( ((beta*u[1]*u[2])/N)^2 +
↪ (((beta*u[1]*u[2])/N) + (gamma*u[2]))^2 + (gamma*u[2])^2 )) du[3] = sqrt(abs( ((beta*u[1]*u[2])/N)*(gamma*u[2]) )) end
At the input, it receives exactly the same parameters as the previous function. () . 4. Obtain a solution: prob = SDEProblem(sir, sir_noise, u0, tspan) sol = solve(prob, EM()) In order to solver . could correctly solve the problem you must specify the type of problem. In this case, for this we use the () function, which receives as input the deterministic form of our model (the () ) function) and its stochastic form (the () function), the initial conditions (0 ) and the time interval ( ) in which diferentiation is performed. Next, we get a solution using the () function, to which the already transformed form of our stochastic model ( ) is fed to the input and a function that implements the solution algorithm is set, in this case we chose the Euler-Maruyama algorithm (function () ).
This completes the solution of the system of stochastic equations. 2.3. Features of solving diferential equations in the Julia programming language In this part of the work, the subtleties that must be taken into account when working with diferential equations in the Julia programming language will be described.
First of all, when working with diferential equations, you need to use the DiferentialEquations.jl package, which contains the necessary tools for solving various kinds of equations. When it comes to obtaining a numerical solution, all the computational work is performed by the () function. For the most eficient construction of work with such a function, general options have been introduced that are compatible with any type of task. In this paper, we will point out some of the useful options: • A standard set of hints for a solution algorithm, or rather a predefined vector _ℎ , which allows you to preset the type of problem.
This is done with the following values: , , . The set value sets the type of the equation as rigid or non-rigid, respectively, the value allows the standard processing algorithm to connect additional procedures to determine the stifness of the problem. The default value is . Additional parameters are available for stochastic equations: ∶ - denotes that the basic SDE has additive noise; ∶ ℎ - Indicates that the solution must be consistent with Stratonovich’s interpretation. • Output control is a set of arguments that control the settings of the solvers, according to which the numerical solution is derived. By default, the maximum set is used to provide the best user experience, but it can be scaled down to save the solution only at a finite point in time.
The parameter designates a specific time to save the numerical solution. The algorithm will store the values at each of the time points in the most eficient way available to the solver. For methods where interpolation is not possible, can be equivalent to .
The parameter denotes the extra time that the time step algorithm must step over to help the solver handle discontinuities and singularities, since going exactly at the break will improve accuracy. If the method cannot change the time steps (multi-step methods with a fixed time step), then the time stops will use interpolation that matches the behavior of . If the method cannot change time steps, and also cannot interpolate, then the value of must be a multiple of , otherwise an error will be thrown. • Step dimension control is a set of arguments controlling procedures associated with setting the time step.
The parameter sets the initial step size. It is also the step size for fixed time step methods. The value is selected automatically if the method is adaptive.
The parameter sets the maximum value of the parameter for the adaptive time step. The default values are package dependent.
The parameter sets the minimum value for the parameter for the adaptive time step. The default values are package dependent.
An additional part in the process of obtaining a numerical solution of equations is the selection of the most suitable algorithm. Some algorithms can give a more accurate solution, while others allow you to get it in a shorter time, but with a more significant error. As an overview, let us designate some of the possible algorithms for application to the types of problems demonstrated in our work.
Speaking about the solving of the ODE system describing the deterministic form of the system, we can consider the following algorithms for its solution: • Euler: Direct Euler method with fixed variable step.; • Midpoint - second order midpoint method. Uses built-in Euler method for adaptability; • Heun: second-order Heun’s method, Euler’s method is used in a similar way; • Ralston - Optimized second-order midpoint method, similarly used Euler’s method; • RK4: Runge–Kutta method of the fourth order. An error control technique is used for an adaptive step with a maximum error over the entire interval; • BS3: Bogatsky-Champin method 3/2.
When reducing our ODE to the form of a jump-like process, we are able to solve such a system using the following algorithm: • SSAStepper: step-by-step algorithm for pure problems like ConstantRateJump and / or MassActionJump. It supports DiscreteCallback handling and persistence control options such as saveat.
To solve the SDE system, we can use the following algorithms: • EM: Euler-Maruyama method. Has a strict order of 0.5 in Ito terms. It can handle all kinds of noise, including of-diagonal, scalar and color noise. The peculiarity of the algorithm is that for its use it is necessary to set only a fixed time step; • EulerHeun: Euler-Heun method. Strong order equal to 0.5 in terms of Stratonovich. It can handle all kinds of noise, including of-diagonal, scalar and color noise. Fixed time step only; • LambaEM: modified Euler-Maruyama method with adaptive time step with error estimation, research authors Lambe and Rakkaukas. Can handle all kinds of noise, including of-diagonal, scalar and color noise; • SRI: the algorithm has a strict order of 1.5 for Ito’s diagonal / scalar SDEs. The backend is based on the SRIW1 method; item SRIW1 - adaptive strong order 1.5 and weak order 2.0 for diagonal / scalar SDEs in Ito form.
Before proceeding to the descriptive part of the materials (graphs) that were obtained in the course of this research work, we want to note that for the deterministic case describing the epidemic model, solution graphs were immediately built without any modification. To obtain graphs describing the behavior of the epidemic model with the addition of discrete stochastic, the kinetic Monte Carlo method (Gillepsy algorithm) built into the DiferentialEquations.jl library was used. And finally, in order to obtain a graph describing the most realistic variant of the model’s behavior (with the addition of self-consistent stochastic), the method of stochastization of one-step processes was used. Thus, in the course of this work, the graphs of the solution for the system of equations of the SIR model were obtained: • In the deterministic form (Fig. 2); • With the addition of discrete stochastic (Fig. 3); • With the addition of self-consistent stochastic (Fig. 4); • General graph comparing all solutions (Fig. 5).
Graph describing the deterministic case solution to the epidemic model (Fig. 2) has a sleek look. This is due to the fact that the system does not have any impacts on the model (external / internal).
Judging by the resulting graph, it can be seen that the number of susceptible individuals (parameter () marked as the blue curve) is rapidly falling and reaches zero by the time of the model time equal to 125 seconds. At the same time, the number of infected individuals (parameter () marked as the red dashed curve) in the same rapid manner begins to increase until the number of susceptible individuals has a value close to zero. This moment comes when the value of the model time approaches 100 seconds. The curve responsible for the designation of the number of afected individuals (parameter () marked as the green dashed curve) increases throughout the model time, dynamically changing the growth rate. The ratio of the growth rate of recovered individuals is associated with the onset of the peak value on the graph of the curve indicating the number of infected individuals (75 seconds of model time). As soon as the number of infected individuals approaches the maximum value and stops growing, the number of infected individuals begins to grow rapidly. Accordingly, when the number of infected individuals approaches zero, the growth dynamics of the curve denoting those who have recovered falls.
Graph describing the solution of the model with the addition of discrete stochastic (Fig. 3) has a similar nature of behavior, however, there are minor fluctuations. This is due to the fact that an impact is added to the system based on a discrete stepwise process.
Graph describing the solution of the stochastized model using the method of stochastization of one-step processes(Fig. 4) has more pronounced fluctuations in comparison with the two previous cases.
This is due to the fact that a term with a self-consistent stochastic is added, which is responsible for external and internal influences. This approach requires clarification: due to the added stochastics, cases of the appearance of new individuals in the system from the vacuum, as well as their disappearance, are possible. Described behavior of model characteristics (() , (), () ) can be seen on the graph at the time interval from 50 to 100 seconds of model time.
There is also a general graph comparing all solutions. (Fig. 5).
In the course of this work, the method of stochastization of one-step processes was applied to the SIR epidemic model.
A software implementation of the computational part for obtaining a solution to the ODE system has been developed. The broad capabilities of the high-level programming language Julia are demonstrated when working with diferential equations of diferent types.
A solution is obtained for each type of model, graphs are constructed that describe the main probabilistic-temporal characteristics of the model, and a numerical analysis of each of the model variants is performed.
This paper has been supported by the RUDN University Strategic Academic Leadership Program and by Russian Foundation for Basic Research (RFBR) according to the research project No 1901-00645.