Neural network approach to data-driven estimation of chemotactic sensitivity in the Keller-Segel model

1.   Introduction

Chemotaxis is the process by which cells or organisms change the state of migration by reacting to the presence of chemicals to access chemically favorable and adverse environments. It plays an important role in many cases, for example, that bacteria move toward the hight concentration of foods and eukaryotic cells migrate towards or away from the chemical source ^[1,2]. Several aspects of chemotactic motility have been studied in great detail, especially for the model organism E. coli. Keller and Segel propsed a mathematical model with chemotaxis in 1970 ^[3]. The model described the behavior of chemotaxis well and has been widely used for the modeling chemotaxis ^[4,5].

This article considers the mathematical model of coinage introduced by Patlak, Keller, and Segel ^[3,6,7,8]. Aggregation and propagation exist everywhere in the population dynamics of chemotactic cells. Chemical chemotactic bacteria are well known for expanding their habitats by patterning them into localized populations. Aggregation occurs in the chemotaxis of a migrating cell that is attracted to move the cell to a higher concentration of the region of the chemical signal it produces. The reaction-diffusion-progression equation for describing the population dynamics of chemotactic cells was proposed by Keller and Segel, where the term of reaction refers to local population dynamics such as cell proliferation and the term of diffusion refers to random movement and convection of cells. A variant of the classic Keller-Segel model was developed to obtain a model that could explain more realistic behavior of chemotactic cells.

Focus on the following equation, the Patlak-Keller-Segel equation.

We study on deep learning architectures to solve Patlak-Keller-Segel equation. Using the initial, boundary condition, and pde equations as loss functions, the Neural Net will be trained as an approximation solution. This approach has two major advantages. First, a function learned in our method can make predictions on any gird. A function trained with a neural net gives an output for any input value. This is a difference from the existing numerical approach that predicted only the value of a specific gird. Second, the unknown coefficient can be predicted by learning. It is sometimes necessary to make predictions when the coefficient is unknown. In this Forward-Inverse Problem, the Neural Net approach is easily applied^[9]. Each coefficients parameterizes as weights that can be learned.

There have been attempts to solve Differential Equation using Neural Nets ^[10,11]. Neural Nets can fit any function ^[12]. In the case of PDEs with complex behaviors, it is difficult to find a solution. While studies to find an approximation solution using a Neural Net in such a PDE has been started recently, there are not many works with the Neural Net approaches in the field of mathematical biology. In this paper, we applied the Neural Net approach to the Patlak-Keller-Segel equation. It is shown in the paper that the Neural Net method can have a good accuracy compared to the existing method, Finite Volume Methods (FVM).

2.   Materials and methods

2.1.   Approximation solution

Keller-Segel equations deal with several numerical approaches. Numerical methods have been implemented using methods such as finite volume methods and finite element methods. We prove existence and uniqueness of a numerical solution to the proposed scheme. Then, we give a priori estimates and establish a threshold on the initial mass, for which we show that the numerical approximation converges to the solution to the PKS system when the initial mass is lower than this threshold. Numerical simulations are performed to verify accuracy and the properties of the scheme. Finally, in the last section we investigate blow-up of the solution for large mass.

In ^[11], a mass transport steepest descent scheme was proposed to solve the modified 1D Keller-Segel system for the log interaction kernel proposed in ^[9]. This method is designed based on the transformation method of the Fokker-Planck type equation introduced in ^[13,14] and applied to the Keller-Segel type model in to satisfy the discrete free energy dissipation principle. By taking into account the problem of transformed variables, this method can accurately solve areas of high concentration without fine-tuning the mesh. This approach has been extended to several dimensions for nonlinear aggregate-diffusion equations and from discretization in ^[15,16] references therein to other approaches.

In these methods, the $\chi$ is entered as a condition. Looking at the CFL condition, it is necessary to hold the grid more closely according to the $\chi$. This is why it is difficult to analyze the Inverse Problem numerically. It is impossible to satisfy the CFL condition with an unknown $\chi$ value. A new method is needed to deal with the inverse problem.

The neural net can be used to find the approximate solution of the PDE. We define the loss function according to equation, boundary and initial condition. Neural net structure is built and learning is performed using the loss function. It has been proved that this approximate solution can always be obtained using the universal approximation theorem. The Neural Net approach is useful for solving inverse problems. We define unknown coefficients as variables.

The neural net approach is useful for solving inverse problems. We define unknown coefficients as variables and proceed with learning. When learning through the optimizer, a coefficient close to the actual value can be obtained. This method does not depend on the value of $\chi$.

Neural network approaches rarely discuss errors between actual and predicted values. It cannot be easily discussed because the actual value is unknown. For this reason, existing studies often measure only the loss used for learning. We compensate for this weakness by bounding the error using loss. The next subsection is about this topic.

2.2.   Deep Learning algorithm

In this section, we introduce our DNN structure and Numerical Method. The proof process in this section was applied by modifying the proof applied to the Fokker-Planck equation.^[17]

Let $f^{NN}(t, x; m, w, b)$ be a deep learning function with smooth activation function. $f$ will be approximation solution of Keller-Segel equation. The structure of f is as follows.

2.3.   Input points

We need the data of grid points for each variable domain. We use random sampling to pick grid points within the domain. The gird points of initial condition and the boundary condition are the same as $\{(t = 0, x_j)\}_j$ and $\{(t_i, x = 1 \text{ or } -1)\}_i$, respectively.

2.4.   Loss functions

We use the Adam optimizer to finds the optimal parameters $w_{ji}^{(l + 1)}$ and $b_j^{(l + 1)}$ to minimize loss functions using the back-propagation method.

The governing equation loss is defined as follows:

Then we define $\text{Loss}_{\text{GE}}$ as

The initial grid points as follows:

The boundary condition is defined as follows:

Finally, we define the total loss as

2.5.   Error bound with Loss

Theorem 1. Let

Then,

for some constant $C$, where $(u^*, v^*)$ is the solution of equation.

Proof. Consider

Integrate on $I$, then

Integrate on $[0, T]$, then

Subtract same equation of $u^*$.

Also consider

Remark 2. The assumption $f \in \widehat{C}^{(1, 1)}([0, T] \times [0, 1])$ can be replaced by a general Sobolev space, since the functions in a Sobolev space can be approximated by the continuous functions on a compact set.

Let $u$ and $v$ be solutions and $u^{\text{NN}}$ and $v^{\text{NN}}$ be approximation solutions. ans $u^d = u - u^{\text{NN}}$, $v^d = v - v^{\text{NN}}$.

Consider

Integrate on $x$, then the first term of RHS

The second term of RHS

Define $B_1 = 2\int_{\partial I} u^d(u\partial_x v)^d n_x dS_x$, and $B_2 = 2\int_I \partial_x u^d (u\partial_x v)^d dx$.

where $||u^d||_{L^\infty(I)} < M_1$.

We may assume $||d_{\text{bd}}^{(2)}||_{L^\infty (\partial I)} < 1/6$.

Choose $\epsilon_1$ such that $1/3 + M_1\epsilon_1 < 1/2$,

Consider $B_2(t)$

First term of RHS is transformed as follows.

where $||u^{\text{NN}}||_{L^\infty(I)} \leq ||u||_{L^\infty(I)} + ||u^d||_{L^\infty(I)} < M_2$. Second term of RHS is transformed as follows.

where $\lVert \partial_x v \rVert_{L^\infty(I)} < M_3$. Hence

Choose $\epsilon_3$ such that $M_2\epsilon_2 + M_3\epsilon_3 < 1/2$,

Also consider,

Silimarly,

Calculate (2)${} + M \times {}$(3),

Hence,

for some constant $C$.

3.   Results

3.1.   Theoretical results

In most cases, $L^2$ loss is used in the PDE solution approximation through Neural Net. As mentioned in the previous subsection, this does not guarantee error bound. We found a loss that can bound the error through the energy estimate.

The energy estimate was calculated from the given PDE equation. By modifying this equation, we proved that the error, the difference between the actual value and the predicted value, is bound to a constant multiple of the loss we are learning. From this proof it can be seen that the approximation with a sufficiently small loss is very close to the actual solution.

The found loss uses $H^1$ loss instead of $L^2$ loss in some terms. This means that reducing the difference in function values is insufficient to reduce errors. If $H^1$ loss is used, not only the function value but also the difference between the derivative value is reduced. In order to make the error sufficiently small, the derivative value of some terms must be small.

Our method is directly applicable to the Inverse Problem. When there is an unknown coefficient, most of the methods cannot be applied as it is. In particular, Numerical methodologies are developed depending on coefficients. It is only possible to apply a new method of predicting parameters through observed data points. On the other hand, our method needs no transformation. If it converges sufficiently, it becomes a good approximator for both the parameter and the solution.

Chemotaxis has been primarily studied as an average characteristic of a population, with little regard for variability among individuals. Despite the simplifying approximations involved in our derivations, especially in the extrapolation to higher spatial dimensions, the models demonstrate a satisfactory and very useful ability to quantitatively interpret population assays for bacterial and leukocyte chemotactic migration.

3.2.   Simulation results

The results from the previous method are as follows. It can be seen that the Neural Net method fits well with the Keller-Segel solution like well-known numerical method. A finite volume method (FVM) was used for the comparison. The simulations were performed using Clawpack ^[18,19] with a fractional step method ^[20] in Python 3.6.

Figure 1 and Figure 2 show the result of Neural Net approximation and FVM approximation, repectively. It can be seen that the Neural Net method is not behind the FVM. Our approximate solution has a chemotaxis phenomenon. One can see that the aggregation takes place over time. This new methodology can sufficiently replace the existing methodology.

We find the coefficients of the solution and the solution itself simultaneously. It couldn't be done with the Numerical method. Find the coefficient and solution simultaneously through the observed data points. It can be seen that there is very little error even when comparing the found coefficient with the actual value.

Our method has the advantage of being able to deal with the inverse problem as well as the forward problem compared to the previous studies. As for the existing numerical method, few studies of the Inverse Problem in Keller-Segel have been conducted. In the case of Neural Net, the forward-inverse problem can be easily solved through simple transformation.

4.   Conclusion

Neural network approaches rarely discuss errors between actual and predicted values. It cannot be easily discussed because the actual value is unknown. For this reason, existing studies often measure only the loss used for learning. We compensate for this weakness by bounding the error using loss.

$H^1$ loss is used instead of $L^2$ loss in some terms. This means that reducing the difference in function values is insufficient to reduce errors. If $H^1$ loss is used, not only the function value but also the difference between the derivative value is reduced. In order to make the error sufficiently small, the derivative value of some terms must be small.

Acknowledgments

Seongwon Lee was supported by National Institute for Mathematical Sciences (NIMS) grant funded by the Korea government (MSIT) (No. B21810000). Hyung Ju Hwang was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. NRF-2017R1E1A1A03070105 and NRF-2019R1A5A1028324), Institute for Information & Communications Technology Promotion (IITP) grant funded by the Korea government (MSIT) (No. 2019-0-01906, Artificial Intelligence Graduate School Program (POSTECH)), and the Information Technology Research Center (ITRC) support program (No. IITP-2018-0-01441).

Conflict of interest

The authors have no conflict of interest.