The classical Lotka-Volterra (L-V) systems are differential equations modeling competitions among multiple species. Research on these systems has been extensive in both theoretical and numerical aspects. Because of uncertain nature of the surrounding environments, L-V systems are subject to random perturbations, and these systems become stochastic as a result. The stochastic L-V system has attracted particular attention in recent years, and studies with interesting results have constantly been emerging. (see Refs.[1-2] and references therein). Related work includes theoretical analysis, numerical simulations, as well as model (parameter) estimations.

With the rapid development of big data science and artificial intelligence^[3], machine learning methods are gradually penetrating into more and more research fields (see Refs. [4-11]), such as the regression of differential equations from observational data using neural networks.

The regression of differential equations can be classified into two categories: deterministic and stochastic. For deterministic ordinary differential equations (ODEs), Raissi et al.^[4] proposed a multi-step network based on the multi-step methods to learn the equations satisfied by the data points. Chen et al.^[5] proposed the neural ODE method, which established the so-called adjoint equations to remove the memory-expensive backpropagating through numerical steps of the ODE integrator. Gholami et al.^[6] gave the ANODE method, where an adjoint based neural ODE framework was designed to avoid problems caused by numerical instability, and unconditionally accurate gradients were provided. Kloppers et al.^[12] discussed the integral and log integral methods for estimating parameters of the L-V models. Wu et al.^[13-14] proposed a linear programming parameter estimation approach according to the grey direct modeling method given. Xu et al.^[15-16] applied the regularization method to inverting parameters, and the inverse problem of parameter estimation from one solution of the system was studied by Shakurov et al.^[17] who also gave the existence and uniqueness conditions. Waniewski et al.^[18] estimated parameters via individual based computer simulations of L-V systems, while Luus^[19] applied direct search optimization method to parameter estimations of L-V systems.

For learning stochastic differential equations (SDEs), Urain et al.^[7] proved the Lyapunov stability of a class of SDEs, and proposed the ImitationFlow which extends the normalizing flows approach to learn stable SDEs from a set of observed trajectories. Kong et al.^[8] gave a neural SDE model (SDE-Net), which consists of a diffusion network and a drift network, learning the diffusion and the drift coefficients respectively. Yildiz et al.^[9] introduced a novel paradigm for learning non-parametric drift and diffusion coefficients of SDEs. The proposed model learns to simulate path distributions that match observations with non-uniform time increments and arbitrary sparseness, which is in contrast with gradient matching that does not optimize simulated responses. Bento et al.^[10] considered linear models for SDEs. We remark that any such model can be associated with a network (i.e. a directed graph) describing which degrees of freedom interact under the dynamics, and tackling the problem of learning such a network from observation of the system trajectory. Tzen et al.^[11] established a connection between infinitely deep residual networks and solutions to SDEs.

The goal of this paper is to estimate learnable parameters of a class of stochastic L-V systems from given data using neural networks. Given observational data from a certain population dynamical system, we first derive the statistics, i.e., the approximate mean and covariance matrix of the population variable based on the Euler-Maruyama discretization scheme. We then design a loss function using the statistics of the data based on the Euler-Maruyama scheme. Finally, we apply the neural network with stochastic gradient descent (SGD) to learn the parameters of the underlying stochastic L-V equations.

1 The neural network algorithm

The underlying stochastic L-V systems that we attempt to learn are of the form^[1]

$\left\{\begin{array}{l} \mathrm{d} \boldsymbol{y}(t)=\boldsymbol{D}(\boldsymbol{y}(t))[(\boldsymbol{b}-\boldsymbol{A} \boldsymbol{y}(t)) \mathrm{d} t+\boldsymbol{\sigma} \mathrm{d} W(t)], \\ \boldsymbol{y}\left(t_0\right)=\boldsymbol{y}_0, \end{array}\right.$

(1)

where y(t)=(y¹, …, yⁿ)^T, and D(y)=diag{y¹, …, yⁿ}) is a diagonal matrix. n represents the number of species in the system, and yⁱ is the population size of the species i. Elements of the vectors b=(b₁, …, b_n)^T and σ=(σ₁, …, σ_n)^T, as well as those of the symmetric matrix A=(a_ij)_n×n are all positive learnable para-meters. W(t) is a 1-dimensional standard Wiener process. For i, j=1, …, n, a_ii denotes the coefficient of intraspecific competition in the species i, a_ij measures the interaction between species i and j, b_i is the intrinsic growth rate of the species i, and σ_i represents the intensity of the noise on the species i.

The Euler-Maruyama scheme^[20] applied to system (1) reads:

$\boldsymbol{y}_1=\boldsymbol{y}_0+\boldsymbol{D}\left(\boldsymbol{y}_0\right)\left[\left(\boldsymbol{b}-\boldsymbol{A} \boldsymbol{y}_0\right) h+\boldsymbol{\sigma} \sqrt{h} \xi\right],$

(2)

where h is the time step size and ξ~N(0, 1).

Denote y(t₀+h) by $\hat{\boldsymbol{y}}_1$. Suppose we are given the observation data {$\left(\boldsymbol{y}_0^{(m)}, \hat{\boldsymbol{y}}_1^{(m, k)}\right)$, m=1, …, N₀, k=1, …, N₁}. Based on the Euler-Maruyama discretization (2), we can approximately regard $\hat{\boldsymbol{y}}_1^{(m, k)}$ as being obtained by

$\begin{gathered} \hat{\boldsymbol{y}}_1^{(m, k)} \approx \boldsymbol{y}_1^{(m, k)}=\boldsymbol{y}_0^{(m)}+ \\ \boldsymbol{D}\left(\boldsymbol{y}_0^{(m)}\right)\left[\left(\boldsymbol{b}-\boldsymbol{A} \boldsymbol{y}_0^{(m)}\right) h+\boldsymbol{\sigma} \sqrt{h} \xi^{(m, k)}\right], \end{gathered}$

(3)

where {ξ^{(m, k)}, m=1, …, N₀, k=1, …, N₁} are independent N(0, 1) random variables.

According to (3), for each m∈{1, …, N₀}, $\hat{\boldsymbol{y}}_1^{(m)}$ approximately obeys the normal distribution with mean $\boldsymbol{\mu}^{(m)}=\boldsymbol{y}_0^{(m)}+\boldsymbol{D}\left(\boldsymbol{y}_0^{(m)}\right)\left(\boldsymbol{b}-\boldsymbol{A} \boldsymbol{y}_0^{(m)}\right) h$ and covariance matrix $\boldsymbol{\Sigma}^{(m)}=h \boldsymbol{D}\left(\boldsymbol{y}_0^m\right) \boldsymbol{\sigma}\left(\boldsymbol{D}\left(\boldsymbol{y}_0^m\right) \boldsymbol{\sigma}\right)^{\mathrm{T}}$. The law of large numbers^[21] implies that, for each m∈{1, …, N₀}the mean and covariance matrix of $\left\{\hat{\boldsymbol{y}}_1^{(m, k)}, k=1, \cdots, N_1\right\}$ approaches $\boldsymbol{\mu}^{(m)}$ and $\mathbf{\Sigma}^{(m)}$ respectively as the number of samples N₁ increases, namely, for sufficiently large N₁

$\begin{array}{l} {{\bf{\mu }}^{(m)}} \approx \overline {\widehat {\bf{y}}_1^{(m)}} = \frac{{\sum\limits_{k = 1}^{{N_1}} {\widehat {\bf{y}}_1^{(m,k)}} }}{{{N_1}}},\\ {{\bf{\Sigma }}^{(m)}} \approx {\bf{\Sigma }}_{{{\hat y}_1}}^{(m)} = \frac{{\sum\limits_{k = 1}^{{N_1}} {\widehat {\bf{y}}_1^{(m,k)}} {{\left( {\widehat {\bf{y}}_1^{(m,k)}} \right)}^{\rm{T}}}}}{{{N_1}}} - \\ \;\;\;\;\;\;\;\;\;\frac{{\sum\limits_{k = 1}^{{N_1}} {\widehat {\bf{y}}_1^{(m,k)}} }}{{{N_1}}}\frac{{{{\left( {\sum\limits_{k = 1}^{{N_1}} {\widehat {\bf{y}}_1^{(m,k)}} } \right)}^{\rm{T}}}}}{{{N_1}}}. \end{array}$

(4)

Thus, for our neural work, we set up the training data $\left\{\left(\boldsymbol{y}_0^{(m)}, \overline{\hat{\boldsymbol{y}}_1^{(m)}}, \boldsymbol{\Sigma}_{\hat{y}_1}^{(m)}\right), m=1, \cdots, N_0\right\}$.

Due to symmetry of the matrix A, there exists a function $\hat{f}(\boldsymbol{y})=\boldsymbol{b}^{\mathrm{T}} \boldsymbol{y}-\frac{1}{2} \boldsymbol{y}^{\mathrm{T}} \boldsymbol{A} \boldsymbol{y}$ such that $\boldsymbol{b}-\boldsymbol{Ay}=\nabla \hat{f}(\boldsymbol{y})$. Meanwhile, it holds that $\sigma=\nabla \hat{g}(\boldsymbol{y})$ for $\hat{g}(\boldsymbol{y})=\sum_{i=1}^n \sigma_i y^i$.

Then we construct the loss function L of our network as follows

$\begin{gathered} \operatorname{Loss}_1=\left\|\overline{\hat{y}_1}-\left(\boldsymbol{y}_0+h \boldsymbol{D}\left(\boldsymbol{y}_0\right) \nabla \hat{f}\left(\boldsymbol{y}_0\right)\right)\right\|_2^2, \\ \operatorname{Loss}_2=\left\|\boldsymbol{\Sigma}_{\hat{y}_1}-h \boldsymbol{D}\left(\boldsymbol{y}_0\right) \nabla \hat{g}\left(\boldsymbol{y}_0\right)\left(\boldsymbol{D}\left(\boldsymbol{y}_0\right) \nabla \hat{g}\left(\boldsymbol{y}_0\right)\right)^{\mathrm{T}}\right\|_2^2, \\ L=\lambda_1 \operatorname{Loss}_1+\lambda_2 \operatorname{Loss}_2\left(\lambda_1, \lambda_2>0, \lambda_1+\lambda_2=1\right), \end{gathered}$

(5)

where the norm ‖·‖₂ denotes the Euclidean norm.

The functions $\hat{f}$ and $\hat{g}$ can then be learned by the net through performing SGD on L. Consequently, we get the estimations

$\begin{gathered} \boldsymbol{A} \approx-\nabla^2 \hat{f}, \\ \boldsymbol{b} \approx \nabla \hat{f}\left(\boldsymbol{y}_0\right)-\nabla^2 \hat{f}\left(\boldsymbol{y}_0\right), \quad \forall \boldsymbol{y}_0, \\ \boldsymbol{\sigma} \approx \nabla \hat{g} . \end{gathered}$

(6)

We write the above process into Algorithm 1.

Algorithm 1: Learning process of stochastic L-V systems
Input: Initial point set $\left\{\boldsymbol{y}_0^{(m)}, m=1, \cdots, N_0\right\}$
Output: The estimate of A, b, σ
1 Getting the observed value {$\hat{\boldsymbol{y}}_1^{(m, k)}, m=1, \cdots, N_0, k=1, \cdots, N_1$};
2 For m=1: N0 do
3   Calculating $\overline{\hat{y}_1^{(m)}}$ and $\Sigma_{g_1}^{(m)}$ according to formula (4);
4 end
5 $\boldsymbol{D} \leftarrow\left\{\left(\boldsymbol{y}_0^{(m)}, \overline{\hat{\boldsymbol{y}}_1^{(m)}}, \Sigma_{\hat{\boldsymbol{y}}_1}^{(m)}\right), m=1, \cdots, N_0\right\}$;
6 Training netNN(·)with data set D;
7 $\hat{f}$←netNN[1];
8 $\hat{g}$←netNN[2];
9 Calculating estimates according to equation (6).

2 Implementation

In this section, we apply our neural network algorithm to the following stochastic L-V system with 2 species^[1], to illustrate the implementation of the algorithm and tesitify its ability of parameter estimation

$\left\{\begin{array}{l} \mathrm{d} x(t)=x(t)\left[\left(b_1-a_{11} x(t)-a_{12} y(t)\right) \mathrm{d} t+\right. \\ \;\;\;\;\left.\sigma_1 \mathrm{~d} W(t)\right], \\ \mathrm{d} y(t)=y(t)\left[\left(b_2-a_{21} x(t)-a_{22} y(t)\right) \mathrm{d} t+\right. \\ \;\;\;\;\left.\sigma_2 \mathrm{~d} W(t)\right], \end{array}\right.$

(7)

Where b₁=1, b₂=1.2, a₁₁=0.3, a₁₂=0.2, a₂₁=0.2, a₂₂=0.4, σ₁=0.4, and σ₂=0.5.

The experiment in this paper is performed in the Windows system, using the PyCharm development environment, and Python 3.7 as the assembly language. In Algorithm 1, we use Pytorch library^[22] to build a feedforward neural network framework net_NN(·) as shown in Fig. 1, where there is one input layer with 2 neurons, one hidden layer with 400 neurons, one output layer with 2 neurons learning $\hat{f}$ and $\hat{g}$ respectively, and the activation functions are all tanh. We set the proportion of the training set to the total data as P_split=0.8, the learning rate as r_learn=0.001, the sample sizes as N₀=1 000, N₁=1 000, and the number of trainings as N_train=2 000.

Hereby we denote (x, y)^T=: y. Given the observation data {$\left(\boldsymbol{y}_0^{(m)}, \hat{\boldsymbol{y}}_1^{(m, k)}\right), m=1, \cdots, N_0, k=1, \cdots, N_0$}, we get the training data {$\left(\boldsymbol{y}_0^{(m)}, \overline{\hat{\boldsymbol{y}}_1^{(m)}}, \right.\left.\boldsymbol{\Sigma}_{y_1}^{(m)}\right), m=1, \cdots, N_0$} by using (4). Then we train the network with the loss function (5) to obtain the outputs $\hat{f}$ and $\hat{g}$. The parameters of the system can be estimated afterwards via (6).

For numerical validation, we generate the observation data in the following way. First we select the data region Ω_d=[0, π]×[0, π] and sample uniformly in it to obtain {y₀^(m), m=1, …, N₀}. Then for each m, we apply the Milstein scheme^[20]

$\begin{gathered} x_{l+1}^{(m, k)}=x_l^{(m, k)}+x_l^{(m, k)}\left(b_1-a_{11} x_l^{(m, k)}-a_{12} y_l^{(m, k)}\right) \delta t+ \\ x_l^{(m, k)} \sigma_2 \xi_l^{(m, k)} \sqrt{\delta t}+\frac{1}{2} x_l^{(m, k)} \sigma_2^2\left(\left(\xi_l^{(m, k)}\right)^2 \delta t-\delta t\right), \\ x_0^{(m, k)}=x_0^{(m)}, \\ y_{l+1}^{(m, k)}=y_l^{(m, k)}+y_l^{(m, k)}\left(b_1-a_{11} x_l^{(m, k)}-a_{12} y_l^{(m, k)}\right) \delta t+ \\ y_l^{(m, k)} \sigma_2 \xi_l^{(m, k)} \sqrt{\delta t}+\frac{1}{2} y_l^{(m, k)} \sigma_2^2\left(\left(\xi_l^{(m, k)}\right)^2 \delta t-\delta t\right), \\ y_0^{(m, k)}=y_0^{(m)}, \end{gathered}$

(8)

with tiny time step δt=0.001 to system (7), with l=0, …, $M-1 \quad\left(M=\frac{h}{\delta t}\right)$ to simulate the observations $\hat{\boldsymbol{y}}_1^{(m, k)} \approx\left(x_M^{(m, k)}, y_M^{(m, k)}\right)^{\mathrm{T}}(m=1, \cdots, \left.N_0, k=1, \cdots, N_1\right)$. Note that, for analyzing statistics of the observation data we have applied the Euler-Maruyama scheme for simplicity, while for simulating these data, we applied here the Milstein scheme for higher accuracy.

In order to observe the impact of different h on the estimation accuracy, we testify on the set of h: h∈{0.01, 0.02, 0.03, 0.05, 0.1, 0.2, 0.4}. Results are given in Section 3.

Due to errors in the training, the estimated values of the parameters may rely on sampling points. To alleviate this issue, $\varOmega_{\mathrm{e}}=\left[\frac{\pi}{2}-\frac{\pi}{2 \sqrt{2}}, \frac{\pi}{2}+\frac{\pi}{2 \sqrt{2}}\right] \times\left[\frac{\pi}{2}-\frac{\pi}{2 \sqrt{2}}, \frac{\pi}{2}+\frac{\pi}{2 \sqrt{2}}\right]$ as subarea is chosen by us from the region Ω_d on which we calculate the averge of the parameters over selected points as the final estimates. The errors between the estimated parameters and the true parameters are measured by

$\varepsilon_p=|p-\hat{p}|,$

(9)

where $\hat{p}$ represents the estimation of the true parameter p.

Note that we take λ₁=λ₂=1/2 in the loss function (5) in our numerical test.

3 Numerical results

Figure 2 illustrates the training and testing loss against the training epoch, for different h, respectively. It can be seen that, both kinds of loss decrease fast as the training epoch approaches near 200. In addition, they also decrease as h decreases, which may be due to the increasing closeness between the true solution and the numerial Euler-Maruyama scheme, based on which the statistics of the data are derived, for smaller h in (3).

Table 1 lists the estimated values of the parameters for h=0.01 with training epoch 2 000, where the errors reach the magnitude 10^-3.

Figure 3 is the lg-lg drawing of the estimation errors of the parameters against h, where the estimation error is defined by (9). The dashed line is of slope 1, which indicates that the estimation errors of all the parameters are approximately of order 1.

In Fig. 4, we compare the x(t) and y(t) trajectories produced by the learned stochastic L-V equations and those by the true SDE system. We see that they coincide well for both x and y. Hereby we take (x₀, y₀)=(1.5, 10), t∈[0, 10], and the time step for numerical evolution is h₀=0.01.

4 Conclusion

We set up a neural network algorithm to learn the unkonwn parameters of a class of stochastic L-V systems. The Euler-Maruyama discretization enables the approximation of the mean and covariance matrix of the observational variable by data statistics, on the basis of which the loss function is established. Numerical tests show the implementation, validate the method, and indicate the possibility of improving learning accuracy by adjusting time step of the numerical scheme. Future research includes learning unknown parameters which are functions of time t, namely b_i(t), a_ij(t), and σ_i(t) (i, j∈{1, …, n}) for non-autonomous stochastic L-V systems.