The Automation of Stochastization Algorithm with Use of SymPy Computer Algebra Library

SymPy computer algebra library is used for automatic generation of ordinary and stochastic systems of differential equations from the schemes of kinetic interaction. Schemes of this type are used not only in chemical kinetics but also in biological, ecological and technical models. This paper describes the automatic generation algorithm with an emphasis on application details.


Introduction
To describe models, where the intensity of the interaction between components depends on the concentration levels of the components, one can use the schemes of kinetic interactions [1][2][3][4]: Here X i represents the quantity of i-th component, the matrices M = [M s i ] and N = [N s i ] are system state matrices which determine the number of interacting components of the system at each stage of the reaction.k + s and k − s are coefficients of interaction in direct and reverse directions, n is system dimension, S is the number of interactions.
Kinetic schemes are often used to model chemical reactions, biological and ecological systems.From these diagrams one can derive a system of ordinary differential equations (ODEs) which gives the evolution of variables X i (t) over time line.
Our research team has generalized the method of stochastisation of deterministic models [5].This method allows to obtain the stochastic model, based on implicit stochastic properties of the system.All necessary information for the stochastisation process can be extracted from the kinetic interaction schemes.For multi-dimensional systems the stochastisation method requires a large amount of routine work.However, it is possible to program all the steps of the method and obtain the drift vector and the diffusion matrix.We carry out the implementation of the algorithm by using the Python language (version 3) and SymPy library for computer algebra.

An algorithm of obtaining ODE and SDE systems
The initial data are the system state matrices M, N and the coefficients of the interaction between system components k + s and k − s .At the first algorithm stage one has to find a matrix R = M T − N T , the columns of which are denoted by vectors r j .Next, one has to calculate the auxiliary values s + j and s − j : .
Using these values one can obtain the drift vector f(x) and the diffusion matrix G(x).The drift vector is calculated by the formula: , and it can be used to write out the system of ODEs To get the system of stochastic differential equations (SDEs) in Itô form, it is necessary to calculate also the diffusion matrix G(x): After that it is possible to obtain the system of Itô SDE: where W is the n-th dimensional Wiener process.

Brief description of the software implementation
For the automation of the above method we used Python [6] with SymPy [7] module.The program takes symbolic vectors X and K and numeric matrices N and M on input, and gives the drift vector and the diffusion matrix on output.For numerical matrices we use NumPy arrays [7].As interface for our module we provide two functions: These functions return the symbolic drift vector and diffusion matrix.It is possible to use built-in SymPy methods to convert f and G to different formats, such as L A T E X formula notation or expression in any supported programming language.SymPy provides support for a large number of programming languages ranging from popular ones such as C/C++, Fortran, Java etc., to new growing languages such as Julia language.

Examples and program verification
Let us consider two examples.The first example is the Verhulst model describing the growth of population in conditions of limited resources (logistic curve), and the second example is the model of brusselator [8].Ilya Prigogine has proposed this model as a simple example of Belousov-Zhabotinsky reaction type.These models are well known and we use them to test our program and to show algorithm in action.

Verhulst model
The model of logistic growth is described by the following scheme: In this case, the matrices N and M turn into vectors: and vectors r turn into scalars r 1 = 1, r 2 = −1.Because the reaction is not reversible, it is enough to compute only s + : Usually the term k 2 x is discarded due to the fact that the x 2 value is large and the x value is small:

The Brusselator model
The Brusselator model is completely determined by the following scheme: Besides the usual reagents X and Y, the inexhaustible reagents A, B, D and E are also involved into reaction.Their number is considered to be time invariant.They cannot be taken into account while constructing the matrices M and N.However, even if we include them in M and N we will obtain correct results.Let use to demonstrate that fact.
Based on the above matrices, our Python functions return the following drift vector and diffusion matrix:

Conclusion
In this paper, we described an algorithm of obtaining ODE and SDE systems from the schemes.Taking into account the inhibitory interactions, will result in algorithm improvement.