Implementation of an analytical-numerical approach to stochastization of one-step processes in the Julia programming language Arseny V. Fedorov1 , Anna O. Masolova1 , Anna V. Korolkova1 and Dmitry S. Kulyabov1,2 1 Peoples’ Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya St, Moscow, 117198, Russian Federation 2 Joint Institute for Nuclear Research, 6 Joliot-Curie, Dubna, Moscow region, 141980, Russian Federation Abstract 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 DifferentialEquations.jl. As a result of this work, the method of stochastization of one-step processes was applied to the SIR epidemic model. Keywords mathematic modeling, differential equations, schemes of kinetics reactions, stochastization of one-step processes, Julia programming language, DSL, Catalyst.jl, DifferentialEquations.jl 1. Introduction The ordinary differential 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 differential 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 Information and Telecommunication Technologies and Mathematical Modeling of High-Tech Systems (ITTMM-2021), Moscow, Russian, April 19–23, 2021 Envelope-Open arsenyvitalievich@yandex.ru (A. V. Fedorov); anna.masolova@gmail.com (A. O. Masolova); korolkova-av@rudn.ru (A. V. Korolkova); kulyabov-ds@rudn.ru (D. S. Kulyabov) GLOBE https://yamadharma.github.io/ (D. S. Kulyabov) Orcid 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) © 2021 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings http://ceur-ws.org ISSN 1613-0073 CEUR Workshop Proceedings (CEUR-WS.org) 92 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 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 DifferentialEquations.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 DifferentialEquations.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 time- dependent process. Equations are one of the ways to represent the dependences of the concen- trations of substances. In this case, the differential 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]. 𝑑𝑆 −𝛽𝑆𝐼 ⎧ = , ⎪ 𝑑𝑡 𝑁 ⎪ 𝑑𝐼 𝛽𝑆𝐼 (1) = − 𝛾𝐼, ⎨ 𝑑𝑡 𝑁 ⎪ ⎪ 𝑑𝑅 ⎩ 𝑑𝑡 = 𝛾 𝐼 , 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 affected individuals at the time 𝑡, 𝑁 = 𝑆(𝑡) + 𝐼 (𝑡) + 𝑅(𝑡), 𝛽 is the coefficient of the intensity of interaction between healthy individuals and sick individuals (which subsequently leads to a disease), 𝛾 is the coefficient of the intensity of recovery of sick individuals. On the example of the SIR epidemic model (1) [4, 5, 6, 7] demonstrate work Catalyst.jl package. This model was chosen due to its simplicity and frequent use in the field of mathematical modeling. It allows for the assessment and comparative analysis of the spread of the disease among individuals. 93 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 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𝐼 , (2) 𝛾 𝐼− → 𝑅. (3) Now, for the model under consideration, we should write down the interaction scheme (2), (3). Then write out the chemical reactions included in the scheme according to the highlighted stages. In the required format for implementation in Catalyst.jl, the diagram is written out as follows: beta, S + I --> 2I gamma, I --> R Using the @reaction_network macro designed to process a set of reactions we set the interaction scheme, it is a complex chemical reaction (a set of chemical reactions): using Catalyst rn = @reaction_network begin beta/N, S + I --> 2I gamma, I --> R end beta N gamma Further when using a bundle of Jupyter Notebook and a graphical illustrator Graphviz we can visualize the recorded interaction scheme: The Catalyst.jl package allows you to generate DE based on a recorded set of chemical reactions. using DiffEqBase, OrdinaryDiffEq u0 = [500.0, 10.0, 1.0] # [S,I,R] at t=0 p = [0.1, sum(u0), 0.01] # [beta, gamma] 94 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 Figure 1: Interaction diagram for the SIR epidemic model tspan = (0.0, 350.0) op = ODEProblem(rn, u0, tspan, p) sol = solve(op, RK4()) # use RK4 ODE solver The Catalyst.jl module has a built-in ability to add stochastic to the model under consideration. First, a sampling algorithm is used, then a jump-like process is generated using the discretized equations, which is then solved. This can be done with the following set of commands: discrete_prob = DiscreteProblem(rn, u0, tspan, p) jump_prob = JumpProblem(rn, discrete_prob, Direct()) sol = solve(jump_prob, SSAStepper()) However, this method does not allow us to obtain an SDE with a self-consistent stochastic part. 2.1. The one-step processes stochastization method In order to obtain the desired stochastic form of the system describing the SIR epidemic model, we will apply the method of stochastization of one-step processes [8, 9, 10]. 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 difference. Such operators are calculated using the following algorithm: 1. Compiled interaction diagram describing the flow of specific reactions in the system; 1 1 2. For the first reaction from the interaction scheme, the operators 𝐼 и 𝐹 will look like this (4): 1 1 𝐼 = (1, 1, 0)𝑇 , 𝐹 = (0, 2, 0)𝑇 . (4) To find these operators, write down the coefficient 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 1 single copy, so we get the operator 𝐼 = (1, 1, 0)𝑇 . Based on the right side of the reaction, 1 we write out the operator 𝐹 where there is a coefficient 2 before the parameter 𝐼, which 1 as a result forms 𝐹 = (0, 2, 0)𝑇 . 95 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 2 2 3. For the second reaction, the operators 𝐼 and 𝐹 (5) are obtained in a similar way: 2 2 𝐼 = (0, 1, 0)𝑇 , 𝐹 = (0, 0, 1)𝑇 . (5) 4. Then we introduce the state change operator 𝑅𝑖 , obtained from the difference between 𝑖 𝑖 1 2 the operators 𝐹 and 𝐼 . For the first and second reactions, the operators 𝑅 , 𝐹 (6) looks like this: 1 2 𝑅 = (−1, 1, 0)𝑇 , 𝑅 = (0, 1, −1)𝑇 . (6) The next step in the implementation of this technique is to calculate the intensities of the transitions + 𝑆𝛼 , − 𝑆𝛼 𝑆 (7), where 𝑖 - is the reaction number, 𝛼 is the number of reactions in the 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 (7) in general form: 𝑛 𝑛 +𝑆 = +𝑘 ∏ 𝜙𝑖! −𝑆 = +𝑘 ∏ 𝜙𝑖! 𝛼 𝛼 𝑖 𝑖𝛼 , 𝛼 𝛼 𝑖 𝑖𝛼 . (7) 𝑖=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 (8): +𝑆 𝑆! 𝐼! 𝑆 ∗ (𝑆 − 1)! 𝐼 ∗ (𝐼 − 1)! =𝛽∗ ∗ =𝛽∗ ∗ = 𝛽𝑆𝐼 . (8) 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 (9) (for this reaction the inverse intensity − 𝑆2 is also does not exist): +𝑆 𝐼! 𝑅! =𝛾∗ ∗ = 𝛾𝐼. (9) 2 (𝐼 − 1)! 𝑅! 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 diffusion matrix 𝐵. General formulas (10) for calculating 𝐴 and 𝐵 are presented below: 𝐴(𝑆, 𝐼 , 𝑅) = + 𝑆𝛼 ∗ 𝑅𝛼 , 𝐵(𝑆, 𝐼 , 𝑅) = + 𝑆𝛼 ∗ 𝑅𝛼 ∗ (𝑅𝛼 )𝑇 . (10) Find 𝐴 and 𝐵 for our model: 2 −𝛽𝐼 𝑆 𝐴 = ∑ + 𝑆𝛼 ∗ 𝑅𝛼 = +𝑆 1 ∗ 𝑅 1 = 𝛽𝐼 𝑆 ∗ (−1, 1, 0)𝑇 + 𝛾 𝐼 ∗ (0, −1, 0)𝑇 = (𝛽𝐼 𝑆 − 𝛾 𝐼) , (11) 𝛼=1 𝛾𝐼 96 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 2 𝐵 = ∑ + 𝑆𝛼 ∗ 𝑅𝛼 ∗ (𝑅𝛼 )𝑇 = + 𝑆1 ∗ 𝑅1 ∗ (𝑅1 )𝑇 + + 𝑆2 ∗ 𝑅2 ∗ (𝑅2 )𝑇 = 𝛼=1 −1 0 𝛽𝐼 𝑆 −𝛽𝐼 𝑆 0 𝛽𝐼 𝑆 ∗ ( 1 ) ∗ (−1, 1, 0)𝑇 + 𝛾 𝐼 ∗ (−1) ∗ (0, −1, 1)𝑇 = (−𝛽𝐼 𝑆 𝛽𝐼 𝑆 + 𝛾 𝐼 −𝛾 𝐼) . (12) 0 1 0 −𝛾 𝐼 𝛾𝐼 Based on the obtained drift vector (11) and the diffusion matrix (12) we obtain an equation in the Langevin form (14): 𝑑𝑥(𝑡) = 𝐴(𝜙)𝑑𝑡 + 𝐵(𝜙)𝑑𝑊 , (13) 𝑆 −𝛽𝐼 𝑆 𝛽𝐼 𝑆 −𝛽𝐼 𝑆 0 𝑑𝑊 1 𝑑 ( 𝐼 ) = (𝛽𝐼 𝑆 − 𝛾 𝐼) 𝑑𝑡 + (−𝛽𝐼 𝑆 𝛽𝐼 𝑆 + 𝛾 𝐼 −𝛾 𝐼) (𝑑𝑊 2 ) . (14) 𝑅 𝛾𝐼 0 −𝛾 𝐼 𝛾𝐼 𝑑𝑊 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. To do this, we will set additional system parameters: beta = 0.1 #disease rate gamma = 0.01 #recovery rate N = sum(u0) Thus, a value for the disease rate 𝛽 was set equal to 0.1. For the recovery rate value 𝛾 value was set to 0.05. Cofficient 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 DifferentialEquations.jl library, the model equations are specified as a function. It receives the following parameters as input: • du - array of technical parameters responsible for differentiation; • u - array of basic parameters described in the section ??; • beta, gamma, N - coefficients described above; 97 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 • 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 differentiation 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 differential equations in the Julia programming language In this part of the work, the subtleties that must be taken into account when working with differential equations in the Julia programming language will be described. First of all, when working with differential equations, you need to use the DifferentialEqua- tions.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 efficient 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 stiffness of the problem. The default value is 𝑎𝑢𝑡𝑜. 98 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 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 efficient 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: 99 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 • 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 off-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 off-diagonal, scalar and color noise. Fixed time step only; • LambaEM: modified Euler-Maruyama method with adaptive time step with error estima- tion, research authors Lambe and Rakkaukas. Can handle all kinds of noise, including off-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. 3. Results 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 DifferentialEquations.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 100 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 Figure 2: Deterministic system solution graph for the SIR model when the value of the model time approaches 100 seconds. The curve responsible for the designation of the number of affected 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). 101 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 Figure 3: System solution graph with the addition of discrete stochastic for the SIR model Figure 4: System solution graph with the addition of self-consistent stochastic for the SIR model 4. Conclusion 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 differential equations of different types. 102 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 Figure 5: General graph of the three solutions for the SIR model 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. 103 Arseny V. Fedorov et al. CEUR Workshop Proceedings 92–104 Acknowledgments 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 19- 01-00645. References [1] J. Bezanson, S. Karpinski, V. B. Shah, A. Edelman, Julia: A fast dynamic language for technical computing (2012) 1–27. arXiv:1209.5145 . [2] D. S. Kulyabov, A. V. Korol’kova, Computer algebra in julia, Programming and Computer Software 47 (2021) 133–138. doi:10.1134/S0361768821020079 . arXiv:2108.12301 . [3] A. V. Fedorov, A. O. Masolova, A. V. Korolkova, D. S. Kulyabov, Application of a numerical- analytical approach in the process of modeling differential equations in the julia language, in: Information Technology, Telecommunications and Control Systems (ITTCS), 2020, volume 1694 of Journal of Physics: Conference Series, IOP Publishing, Innopolis, 2020, pp. 012026.1–8. doi:10.1088/1742- 6596/1694/1/012026 . [4] R. Ross, An application of the theory of probabilities to the study of a priori pathome- try.—part i, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 92 (1916) 204–230. doi:10.1098/rspa.1916.0007 . [5] R. Ross, An application of the theory of probabilities to the study of a priori pathome- try.—part ii, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 93 (1917) 212–225. doi:10.1098/rspa.1917.0014 . [6] R. Ross, An application of the theory of probabilities to the study of a priori pathome- try.—part iii, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 93 (1917) 225–240. doi:10.1098/rspa.1917.0015 . [7] W. O. Kermack, A. G. McKendrick, A contribution to the mathematical theory of epidemics, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 115 (1927) 700–721. doi:10.1098/rspa.1927.0118 . [8] A. V. Korolkova, D. S. Kulyabov, One-step stochastization methods for open systems, EPJ Web of Conferences 226 (2020) 02014.1–4. doi:10.1051/epjconf/202022602014 . [9] M. Hnatič, E. G. Eferina, A. V. Korolkova, D. S. Kulyabov, L. A. Sevastyanov, Operator approach to the master equation for the one-step process, EPJ Web of Conferences 108 (2016) 02027. doi:10.1051/epjconf/201610802027 . arXiv:1603.02205 . [10] E. G. Eferina, M. Hnatich, A. V. Korolkova, D. S. Kulyabov, L. A. Sevastianov, T. R. Velieva, Diagram representation for the stochastization of single-step processes, in: V. M. Vish- nevskiy, K. E. Samouylov, D. V. Kozyrev (Eds.), Distributed Computer and Communication Networks. DCCN 2016. Communications in Computer and Information Science, volume 678 of Communications in Computer and Information Science, Springer International Publishing, Cham, 2016, pp. 483–497. doi:10.1007/978- 3- 319- 51917- 3_42 . 104