Linearized Implicit Methods Based on a Single-Layer Neural Network: Application to Keller-Segel Models

04/08/2020 ∙ by M. Benzakour Amine, et al. ∙ 0

This paper is concerned with numerical approximation of some two-dimensional Keller-Segel chemotaxis models, especially those generating pattern formations. The numerical resolution of such nonlinear parabolic-parabolic or parabolic-elliptic systems of partial differential equations consumes a significant computational time when solved with fully implicit schemes. Standard linearized semi-implicit schemes, however, require reasonable computational time, but suffer from lack of accuracy. In this work, two methods based on a single-layer neural network are developed to build linearized implicit schemes: a basic one called the each step training linearized implicit (ESTLI) method and a more efficient one, the selected steps training linearized implicit (SSTLI) method. The proposed schemes make use also of a spatial finite volume method with a hybrid difference scheme approximation for convection-diffusion fluxes. Several numerical tests are performed to illustrate the accuracy, efficiency and robustness of the proposed methods. Generalization of the developed methods to other nonlinear partial differential equations is straightforward.



There are no comments yet.


page 12

page 14

page 16

page 17

page 18

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

One of the fundamental characteristics of living organisms is their ability to move in response to external signals. Chemotaxis, which is the oriented motion induced by chemical gradients, has a critical role in many biological and medical fields. It is encountered in several self-organization phenomena such as: aggregation of bacteria, embryogenesis, ecology, tumor growth, etc.

In 1970, Keller and Segel Keller1 introduced a chemotaxis model to describe the aggregation of slim molds. Since then, several variants of the Keller-Segel model have been developed to model phenomena in a variety of fields where chemotaxis is involved (see for more details Painter ). The model is still very popular, it reads in a generalized form


where denotes the cell density at spatial position and time , is the chemoattractant concentration, and respectively define cells diffusion and chemical diffusion coefficients, is the chemotactic sensitivity, is a function describing cell growth and death, and is a kinetic function that describes production and degradation of the chemical signal.

Several numerical methods for Keller-Segel systems have been proposed and analyzed. The interested reader may be referred for instance to the introduction of Akhmouch1 , and to Jungel ; Huang ; Xiao ; Zhang1 ; Sulman ; Zhang2 ; Liu ; Li ; Dehghan for more recent works.

In this paper, we focus on three models derived from (1.1). The first one has been proposed in the framework of pattern formation in embryology, the model reads:


with zero-flux boundary and initial conditions


Here and in the rest of this paper, and are positive constants, , is an open bounded subset, denotes the final time, is the outward unit-normal on the boundary , and is a nonnegative function.

The above model is a simpler version of the original parabolic-parabolic model proposed by Oster and Murray Oster1 , in which the second equation of the system is elliptic, using the reasonable assumption that the chemoattractant diffuses much faster than cells. In Murray1 , it has been demonstrated that the model is able to generate stripe patterns similar to that of embryonic American alligators.

In Akhmouch2 , a semi-implicit finite volume scheme with an additional correction term has been proposed by Akhmouch and the author to solve the system (1.2)–(1.3). Moreover, a convergence analysis has been performed in the case of a nonnegative initial cell density , and it has been established that the numerical solution of the proposed scheme converges to a weak solution of (1.2)–(1.3) under some conditions. In this paper, a particular attention will be given to this analysis.

The second Keller-Segel system considered in this paper is the following chemotaxis-growth model:


with zero-flux boundary conditions and suitable initial conditions. Herein, is a positive constant and is a polynomial growth term. This model has the ability to reproduce different bacterial pattern formations through the coupled effect of chemotaxis and growth (see, e.g., Aida ; Strehl1 ; Strehl2 ; Akhmouch1 ; Akhmouch3 ). It is well known that when , the solution of the system (called minimal model) may blow-up in finite time for initial cell density with sufficiently large mass, which is an interesting mathematical phenomenon. Several works were devoted to study this system (see for more details Horstmann1 ; Hillen1 ; Perthame1 ). However, from biological point of view, the occurrence of finite-time blow-up is considered as a pathological behavior of the model. An adequate choice of the reaction term can prevent this blow-up Winkler .

The last Keller-Segel system considered in the current paper reads:


endowed with zero-flux boundary conditions and initial conditions. In the above model, is a function which satisfies the conditions: , there exists such that and for . The system (1.5), called volume-filling chemotaxis model, has been introduced by Hillen and Painter in Hillen2 , where the existence of solutions that are bounded is proved. The adopted choice of chemotactic sensitivity implies that cells stop to accumulate when the threshold is reached by , which prevents overcrowding. Numerical simulations of the model are performed in Hillen3 , and in Andreianov1 ; Chamoun1 for more general diffusion term.

The three chemotaxis models presented are highly nonlinear coupled systems due to the chemotaxis term ( for the volume-filling model). Moreover, the cell density reaction term of the second equation of (1.2), the growth term of the first equation of (1.4), and the convective term of the first equation of (1.5), are also nonlinear. These nonlinearities are an important factor for the choice of the time discretization approach when numerically solving such systems. The most used approaches in the literature are: explicit, implicit and semi-implicit schemes. The main drawback of explicit schemes is the too restrictive CFL condition necessary to ensure stability of the solution. For implicit schemes, the CFL constraint is avoided. But, on the other hand, a large nonlinear algebra system must be solved at each time-step, which is too time-consuming. The third time discretization which is also widely used for solving nonlinear systems is the semi-implicit discretization. Semi-implicit schemes are linearized implicit schemes, so they only require to solve decoupled linear algebra systems at each time level. Regarding numerical methods for Keller-Segel systems, Euler semi-implicit scheme is used in several works to obtain fully linear implicit schemes (see, e.g., Strehl2 ; Akhmouch3 ; Liu ; Saito1 ; Saito2 ). In other works, it has only been used to avoid the severe nonlinearity caused by the chemotaxis term and thus obtain decoupled nonlinear schemes Andreianov1 ; Chamoun1 ; Ibrahim1 . As proved in the convergence analysis carried out by Zhou and Saito Saito2 for the parabolic-elliptic minimal model, by Akhmouch and the author Akhmouch3 for parabolic-elliptic chemotaxis growth model and by Andreianov et al. for degenerate volume-filling chemotaxis model Andreianov1 , no CFL condition is required to obtain the convergence of Euler semi-implicit finite volume schemes for such systems. Nonetheless, despite the advantages of such strategy, these schemes generally suffer from a lack of accuracy.

In this work, a novel approach to develop linear implicit schemes is presented. More precisely, two methods relying on a basic single-layer neural network to linearize schemes are proposed: the each step training linearized implicit (ESTLI) method, and the selected steps training linearized implicit (SSTLI) method which is based also on a selected steps training (STP) algorithm. The developed methods will be used to numerically solve the systems (1.2), (1.4) and (1.5). The objective is to build a method which has the following attractive benefits:

  • It must only require to solve decoupled linear systems at each time-step.

  • It should be significantly more accurate than semi-implicit method.

  • In terms of computational cost, the method should be nearly as efficient as semi-implicit method, specially for time-consuming numerical experiments.

  • As semi-implicit method, it is required to be stable even for large time-step sizes.

  • It should be easily applied to any equation or system which can be solved by semi-implicit scheme.

Concerning the spatial approximation, a finite volume scheme similar to that of Akhmouch and Benzakour Akhmouch2 is adopted. We mention that finite volume method has been first introduced by Filbet Filbet1 for Keller-Segel models, where a fully implicit finite volume scheme has been proposed for the parabolic-elliptic minimal model.

As mentioned above, the approach presented in this article uses artificial neural networks (ANN) only as a linearization technique in a numerical scheme. We mention, however, that several interesting works were devoted to use ANN as a main tool to approximate solutions of partial differential equations, so we briefly review some of them. In Lagaris1

, multi-layer perceptron was used by Lagaris et al. to solve boundary value problems. Lagaris et al.


developed also a method based on multi-layer perceptron and radial basis function networks. In

Shirvany , a novel method was presented using the same networks, with application to the nonlinear Schrodinger equation. We note that several works proposed radial basis function networks to solve elliptic equations, see for instance Mai-Duy ; Aminataei ; Jianyu . The reliability of other methods based on single-layer Bernstein neural network Sun and single-layer Chebyshev neural network Mall , was also demonstrated by application to elliptic equations. In Rudd , a method combining Galerkin methods and ANN was developed for the approximation of the solution of hyperbolic and parabolic partial differential equations. ANN were also used to solve systems of partial differential equations Beidokhti , and high dimensional partial differential equations Sirignano ; Weinan .

This paper is organized as follows. In Sect. 2, a semi-implicit Euler scheme in time and finite volume in space for the the initial-boundary value problem (1.2)–(1.3) is presented. The convergence of the scheme is also discussed in this section on the basis of the analysis performed in Akhmouch2 . Sect. 3 introduces ESTLI and SSTLI methods to solve (1.2)–(1.3). The STP algorithm is also presented in this section. Numerical tests applied to the models presented in this section are performed in Sect. 4, in order to compare the efficiency and accuracy of the elaborate schemes with that of semi-implicit schemes. Finally, Sect. 5 presents a conclusion.

2 The semi-implicit finite volume scheme

2.1 Definitions and notations

We assume that is an open bounded polygonal subset. An admissible finite volume mesh of in the sens of Definition 9.1 in Eymard1 is given by:

  • A family of control volumes (disjoint open and convex polygons) denoted by .

  • A family of edges, where is the set of edges of the control volume .

  • A family of points such that . The particularity of an admissible mesh is that the straight line going through and must be orthogonal to the common edge of and denoted .

We denote by m, the Lebesgue measure in or . Denoting by d the Euclidean distance, and for all , we define by:

Let be the number of time-steps to reach the final time . We set , where is the time-step size: .

Let be the set of functions from to which are constant over each control volume of the mesh. Let , for , the classical discrete norm reads

where for all and for all . We define also the discrete seminorm and the discrete norm:

where for all , if and otherwise, with .

Finally, we define a weak solution of the system (1.2)–(1.3):

Definition 2.1.

A weak solution of (1.2)–(1.3) is a pair of functions which satisfy the following identities for all test functions :


2.2 Presentation of the scheme

A semi-implicit finite volume scheme associated to the problem (1.2)–(1.3) is given by:

for all and ,


endowed with the discrete initial condition


The function is defined by


where is a small nonnegative constant (), and if and otherwise.

In the above scheme, and denote respectively the approximations of the mean value of and on at time . As we can see, the difference between this scheme and a fully implicit one is the discretization of the reaction term , which is evaluated at the previous time-step. This strategy allows us to build a decoupled linear scheme: at each time-step, we solve (2.4) to obtain and then, we solve (2.3) to compute . The discretization used for the chemotaxis term is identical to the first order upwind scheme when , and to the second order central difference scheme otherwise. It is equivalent to that of Spalding Spalding1 when .

2.3 Convergence of the scheme

The scheme (2.3)–(2.5) is similar to that of Akhmouch2 , where a corrected decoupled scheme is proposed for the problem (1.2)–(1.3). The only difference between the two schemes is a bounded term which was added to the second equation of the scheme:

we refer to Section 2.2 of Akhmouch2 for the definition of . The numerical analysis of the two schemes is quite similar, so based on the results obtained for the corrected decoupled scheme, we briefly discuss the convergence of (2.3)–(2.5).

The existence, uniqueness and nonnegativity of the finite volume solution of (2.3)–(2.5) can be proved by an M-matrix analysis, following exactly the proof of Proposition 3.1 in Akhmouch2 .

Now, to prove the convergence of the scheme to a weak solution of (1.2)–(1.3

) in the sense of Definition 2.1, we first need a priori estimates on the discrete solutions. For the finite volume approximation of

, and discrete estimates can be obtained directly from the equation (2.3). We refer to Proposition 4.1 of Akhmouch1 , where the same equation is analyzed.

For the cell density finite volume solution, only estimate is obtained for the corrected decoupled scheme without establishing the boundedness of the solution, since it is supposed that . Two approaches have been proposed to obtain this estimate. The first one, detailed in the proof of Lemma 4.1 in Akhmouch2 , needs that , and holds also for the semi-implicit scheme (2.3)–(2.5). The main tool of this approach is the discrete Gagliardo-Nirenberg-Sobolev inequality Chatard2 , which requires the following constraint on the mesh: there exists such that The second approach (Proposition 4.3 in Akhmouch2 ) is based on the discrete Gronwall inequality, and requires a time-step condition which can be relaxed in the case of the current semi-implicit scheme.

In the following proposition, estimate on the approximate solution of the semi-implicit scheme is established if is a bounded function, which is the most relevant case for numerical purposes.

Proposition 2.0.

Assume that there exists such that . Assume that is a bounded function. Let be the solution of the scheme (2.3)–(2.5). Then there exists a constant , only depending on and , such that for all , and


Let be the control volume for which , which implies that for all . Using the identity , we can easily see that for all , we have

From (2.3), and using the above equality, we find that


From (2.6), we can easily see that if . In the case when , we have

In view of the choice of , we deduce that

Next, multiplying (2.4) by , we have

Collecting the previous inequalities we obtain

Using the facts that is nonnegative, and is a bounded nonnegative function, the desired result follows by induction on .

The estimate for the cell density can be now obtained by following the proof of Proposition 4.2 in Akhmouch2 . The rest of the convergence analysis of the scheme (2.3)–(2.5), including the compactness of the sequence of approximate solutions and the pass to the limit in the scheme, can be performed exactly as in Section 5 in Akhmouch2 (see also Hillairet1 ; Filbet1 ).

3 Linearization via a single-layer neural network

In this section, for all family of elements , we denote by

the vector whose components are the elements of the family.

3.1 The ESTLI method

The idea of the ESTLI method is to replace (2.4) in the semi-implicit scheme by the equation:


where is an approximation of , obtained through the use of an ANN. The training data are obtained from the previous values of the finite volume solution, which implies that the numerical solution must be computed by an other method for the first time-steps.

It is easy to see that the ESTLI scheme consumes more time than the semi-implicit scheme (2.3)–(2.4). Indeed, at each time step, in addition to the time required for forming and solving linear systems, a computational time is needed to train the neural network and to compute for all

. As training single-layer neural networks already requires much computational time, and keeping in mind that the scheme must be efficient (see the objectives listed in Sect. 1), hidden layers are avoided. The network proposed is a single-layer feedforward neural network which consists of two inputs, weights, a bias and an output. The activation function will be the identity activation function.

The structures of the ANN used in the training phase and prediction phase are presented in Figs. 1 and 2 respectively. As we can see, at each time-step with , the network is trained using the two inputs and

, and the target output

. The semi-implicit scheme (2.3)–(2.4) will be used for the two first time-steps ( and ) to obtain the numerical solution. The values of the discrete solution at each control volume are used in the training process. After completing the training, the final weights and the bias are used to compute , so for all :

Figure 1: Structure of the neural network at time-step , (Training phase)
Figure 2: Structure of the neural network at time-step , (Prediction phase)

The training process adopted is a standard one: the weights and the bias are updated until the error reaches the desired threshold . Denoting by the vector of weights and the bias at the time step and the iteration (related to the updating process), the update rule has the following form:


where the expression of is determined according to the method selected as learning algorithm: gradient descent method, conjugate gradient method, Levenberg-Marquardt method, etc. We refer to Yu ; Tan for an overview of update rules related to several optimization methods. The error function adopted to evaluate the training process is the mean squared error function :


where is the computed output at the time step and the iteration .

One of the advantages of using ANN to build implicit linear schemes is the compatibility between this approach and the implicit step-by-step resolution. Indeed, if for we need to generate a randomly initial vector of weights and the bias, we will set starting from the fourth time-step. Here, is the vector of weights and the bias at the time-steps . This is because it is expected that the behavior of the solution does not change much between two successive time-steps, and therefore will be very close to . This will reduce the number of iterations needed to obtain the final bias and weights, specially for small time-steps.

3.2 The SSTLI method

The SSTLI method is an efficient version of the ESTLI method in which training of the network does not occur at each time-step, but on selected time-steps. The SST algorithm is developed to this purpose. Two criteria must be satisfied to train the network at the time step :

  1. The performance required is not reached for the time-step : It is widely expected that the approximation of obtained by the use of ANN is better than simply use of a semi-implicit linearization, which means that ( is defined by (3.2) and is computed by the scheme (2.3)–(3.1)). To measure the performance of the scheme, the approach proposed is to introduce a parameter , and depending on the choice of this parameter, the scheme (2.3)–(3.1) will be said to satisfy the required performance at the time-step if:


    In this case, the network is not trained for the following time-step, and we set in order to compute .

  2. is a better approximation for than : Defining by, for all , , the condition must be satisfied to train the network at the time-step . Otherwise, will be used to compute and .

The detailed steps of the SST algorithm are presented in the following. The algorithm stops when the the number of time-steps to reach the final time is exceeded.

  1. For and , compute the numerical solution using the semi-implicit scheme (2.3)–(2.4).

  2. For , randomly generate , train the neural network (Fig. 1) until the error , use the obtained vector of weights and bias to compute , and then solve the scheme (2.3)–(3.1) to have . Finally, set and go to the step 4.

  3. If , set and go to the step 6. Otherwise, go to the next step.

  4. If , set and go to the step 9. Otherwise, set , train the neural network (Fig. 1) until the error and obtain , go to the next step.

  5. Compute and defined by: for all , . Solve the scheme (2.3)–(3.1) to obtain , set and go to the step 3.

  6. Compute , solve the scheme (2.3)–(3.1) to obtain , set and go to the next step.

  7. Set . If go to the step 9. Otherwise go to the next step.

  8. Compute , solve the scheme (2.3)–(3.1) to obtain , and set . Then, set , train the neural network (Fig. 1) until the error and obtain a new . Compute a new , set and go to the step 3.

  9. Compute , solve the scheme (2.3)–(3.1) to obtain . Set and go to the step 4.

4 Numerical experiments

This section mainly deals with the comparison between the Euler semi-implicit method, the ESTLI method, and the SSTLI method. In all numerical tests presented in this section, target outputs are normalized to the range of ( by considering the positivity of the numerical solution). The Levenberg-Marquard method is chosen as learning algorithm, and the threshold mean squared error required is . Finally, for the parameter defined in (3.5), we take .

4.1 Keller-Segel model for embryonic pattern formation

In this test, we consider the system (1.2)–(1.3). The spatial domain is with a uniform mesh of control volumes, and the final time is . We adopt the following initial condition:



is a positive perturbation which is constant on each control volume. It is equal on each cell of the mesh to the average of ten uniformly distributed random values in

. For the model parameters, we take: and .

The numerical cell density obtained by the SSTLI method with is presented in Fig. 3. We can observe from this figure the ability of the studied model to generate stripe patterns. We can see also from the three-dimensional plot on the right that the obtained numerical solution is positive and is oscillation free. In Fig. 4, we present three-dimensional plots of the cell density using SSTLI and ESTLI methods with a very large time-step size . As we can see, both numerical solutions are free from negative values or any instabilities.

Figure 3: Solution of (1.2)–(1.3) at final time computed via the SSTLI method with
Figure 4: Solution of (1.2)–(1.3) at final time computed via the SSTLI method (left) and the ESTLI method (right) with

In order to compare the accuracy of the three approaches presented in this paper, and since the exact solution of the Keller-Segel system (1.2)–(1.3) is unavailable, we use a reference solution computed by the corrected decoupled scheme Akhmouch2 on a very fine time-stepping . This reference solution is used to compute the relative -errors at final time (see Table 1) for several time-step sizes. Computational cost (CPU) of the semi-implicit scheme (2.3)–(2.5), and the increase in computational cost (denoted as ) of SSTLI and ESTLI methods over the semi-implicit scheme are also presented in Table 1. From this table, we can observe that SSTLI method is slightly more accurate than ESTLI method except for . Moreover, both methods are in general four to five times more accurate than Euler semi-implicit scheme. Concerning the computational time, the presented results show that ESTLI method is considerably more time-consuming than semi-implicit scheme. On the contrary, the increase in computational cost over semi-implicit scheme is generally weak in the case of SSTLI method, specially for the three finest time-steps (less than ). Finally, in Fig. 5, it can be seen that the three numerical approaches are first-order accurate.

-error () -error ESTLI () -error CPU (s)
SSTLI ESTLI semi-implicit
Table 1: Computational time and relative -errors obtained for the solution of (1.2)–(1.3) using SSTLI, ESTLI, and Euler semi-implicit methods
Figure 5: Convergence speed of the SSTLI, ESTLI, and Euler semi-implicit methods applied to the the model (1.2)–(1.3)

4.2 Keller-Segel model with quadratic growth

In this subsection, we are concerned with the system (1.4) endowed with zero-flux boundary conditions. Herein, the computational domain is with a uniform mesh grid ., , , , and . The following initial conditions are considered:

and . The perturbation is defined as in (4.1). The scheme adopted is:


when using the semi-implicit method, is replaced by , for ESTLI method is computed by an ANN similar to that of Figs. 1 and 2. In the case of SSTLI method, SST algorithm is used.

The numerical cell density computed via SSTLI method with is plotted in Fig. 6. The continuous rings observed in the figure are similar to the patterns formed by Salmonella typhimurium bacteria Woodward1 . Numerical solutions computed by the SSTLI and ESTLI methods with are shown in Fig. 7. We can see that even for a large time-step, the stability and the positivity of the solution are assured.

Figure 6: Solution of (1.4) with quadratic growth at final time computed via the SSTLI method with
Figure 7: Solution of (1.4) with quadratic growth at final time computed via the SSTLI method (left) and the ESTLI method (right) with

The reference solution used is computed using the corrected decoupled scheme with . For a sequence of decreasing time stepping, relative -errors and computational costs of the studied numerical methods are presented in Table 2, convergence speeds are presented in Fig. 8. Table 2 shows that SSTLI and ESTLI methods are much more accurate than semi-implicit method. The obtained results do not allow us to conclude which of the two proposed methods is more accurate, but it can be observed from the table that SSTLI give better results for the two smallest time-steps. In terms of computational cost, the increase in computational time over semi-implicit scheme is very marginal for the three most time-consuming tests when using SSTLI method (less than ), which is not the case for the ESTLI method.

-error () -error () -error CPU (s)
SSTLI ESTLI semi-implicit
Table 2: Computational time and relative -errors obtained for the solution of (1.4) with quadratic growth using SSTLI, ESTLI, and Euler semi-implicit methods
Figure 8: Convergence speed of the SSTLI, ESTLI, and Euler semi-implicit methods applied to the the model (1.4) with quadratic growth

4.3 Keller-Segel model with volume-filling

Here, we focus on the model (1.5) with zero-flux boundary conditions. The computational domain and the initial chemoattractant concentration are taken as in the previous subsection, and we prescribe the following initial cell density

is defined as in (4.1). Besides, we set and . The following scheme is used:

where the value of is determined as in (4.2)–(4.3), and when .

Numerical solution computed by SSTLI method with is plotted in Fig. 9. We observe from the figure that the maximal density of cells does not exceed even if the solution concentrates at the center, which is consistent with our choice of the function .

Figure 9: Solution of (1.5) at final time computed via the SSTLI method with

To compare relative errors for the three numerical methods studied in this paper, the reference solution of this test is computed using the SSTLI method with . Table 3 shows that relative -errors for SSTLI and ESTLI methods are very close with a slight advantage for SSTLI method in the case of small time-steps. The two methods are about three to four times more accurate than semi-implicit scheme, and the three numerical approaches are first-order accurate (see Fig. 10). Table 3 also shows that, except for , the increase in computational cost over semi-implicit scheme is weak to insignificant for SSTLI method. In terms of -error, we observe in Table 4 that SSTLI and ESTLI methods are four to six times more accurate than semi-implicit scheme.

-error () -error () -error CPU (s)
SSTLI ESTLI semi-implicit
Table 3: Computational time and relative -errors obtained for the solution of (1.5) using SSTLI, ESTLI, and Euler semi-implicit methods
Figure 10: Convergence speed of the SSTLI, ESTLI, and Euler semi-implicit methods applied to the the model (1.5)
-error -error -error
SSTLI ESTLI semi-implicit
Table 4: Relative -errors obtained for the solution of (1.5) using SSTLI, ESTLI, and Euler semi-implicit methods

4.4 Keller-Segel model with cubic growth

The purpose of this test is to investigate the ability of SSTLI and ESTLI methods to capture some bacterial patterns which can be reproduced by the model (1.4). More precisely, we focus on bacterial honeycomb pattern, reported for example in Thar , and symmetrical spot pattern generated for instance by Escherichia coli budrene2 .

In this subsection, , and . Initial conditions and computational domain are the same as in Subsection 4.2. The scheme adopted is similar to (4.2)–(4.3), but since the logistic source term is now cubic, in (4.2) will be replaced by . The time-step size used for computations is .

Figure 11: Solution of (1.4) with cubic growth at computed via the SSTLI method (left) and the ESTLI method (right) with (black indicates low cell density, white indicates high cell density)
Figure 12: Three-dimensional plot of the solution of (1.4) with cubic growth at computed via the SSTLI method (left) and the ESTLI method (right) with

For , the computed solutions at are shown in Fig. 11. The honeycomb pattern is observed for both ESTLI and SSTLI numerical solutions. When , we see from Fig. 12 that the computed solutions exhibit a very spiky behavior at , which means that spot pattern appears (see Fig. 13). Despite this fact, the two methods remain positivity preserving which confirms their robustness and reliability.

Figure 13: Solution of (1.4) with cubic growth at computed via the SSTLI method (left) and the ESTLI method (right) with (black indicates low cell density, white indicates high cell density)

5 Conclusion

In this paper, two linearized implicit methods relying on a single-layer neural network are developed to solve two-dimensional Keller-Segel systems: ESTLI and SSTLI methods. In terms of accuracy, the numerical tests performed demonstrate the superiority of the developed methods on semi-implicit method. However, concerning computational cost, only SSTLI method is nearly as efficient as semi-implicit scheme for time-consuming numerical tests. Moreover, the two numerical methods are robust and easily applicable to any problem which can be solved by the semi-implicit approach.


  • (1) Keller, E.F., Segel, L.A.: Initiation of slime mold aggregation viewed as an instability. J. Theor. Biol. 26, 399–415 (1970)
  • (2) Painter, K. J.: Mathematical models for chemotaxis and their applications in self-organisation phenomena. J. Theor. Biol. 481, 162–182 (2019)
  • (3) Akhmouch, M., Benzakour Amine, M.: A time semi-exponentially fitted scheme for chemotaxis-growth models. Calcolo 54(2), 609–641 (2017)
  • (4) Jüngel, A., Leingang, O.: Blow-up of solutions to semi-discrete parabolic-elliptic Keller-Segel models. Discrete Contin. Dyn. Syst. Ser. B, 609–641 (2019)
  • (5) Huang, X., Xiao, X., Zhao, J., Feng, X.: An efficient operator-splitting FEM-FCT algorithm for 3D chemotaxis models. Eng. Comput. 1–12 (2019)
  • (6) Xiao, X., Feng, X., He, Y.: Numerical simulations for the chemotaxis models on surfaces via a novel characteristic finite element method. Comput. Math. Appl. 78(1), 20–34 (2019)
  • (7) Zhang, Y., Zhang, J.: The splitting mixed element method for parabolic equation and its application in chemotaxis model. Appl. Math. Comput. 313, 287–300 (2017)
  • (8) Sulman, M., Nguyen, T.: A Positivity Preserving Moving Mesh Finite Element Method for the Keller-Segel Chemotaxis Model. J. Sci. Comput. 80(1), 649–666 (2019)
  • (9) Zhang, R., Zhu, J., Loula, A. F., Yu, X.: Operator splitting combined with positivity-preserving discontinuous Galerkin method for the chemotaxis model. J. Comput. Appl. Math. 302, 312–326 (2016)
  • (10) Liu, J. G., Wang, L., Zhou, Z.: Positivity-preserving and asymptotic preserving method for 2D Keller-Segel equations. Math. Comp. 87(311), 1165–1189 (2018)
  • (11) Li, X. H., Shu, C. W., Yang, Y.: Local discontinuous Galerkin method for the Keller-Segel chemotaxis model. J. Sci. Comput. 73, 943–967 (2017)
  • (12) Dehghan, M., Abbaszadeh, M., Mohebbi, A.: A meshless technique based on the local radial basis functions collocation method for solving parabolic-parabolic Patlak-Keller-Segel chemotaxis model. Eng. Anal. Bound. Elem. 56, 129–144 (2015)
  • (13) Oster, G.F., Murray, J.D.: Pattern Formation Models and Developmental Constraints. J. expl. Zool. 251, 186–202 (1989)
  • (14) Murray, J. D., Deeming, D. C., Ferguson, M. W. J.: Size-dependent pigmentation–pattern formation in embryos of Alligator mississippiensis: time of initiation of pattern generation mechanism. Proc. R. Soc. B 239, 279–293 (1990)
  • (15) Akhmouch, M., Benzakour Amine, M.: A corrected decoupled scheme for chemotaxis models. J. Comput. Appl. Math. 323, 36–52 (2017)
  • (16) Aida, M., Tsujikawa, T., Efendiev, M., Yagi, A., Mimura, M.: Lower estimate of the attractor dimension for a chemotaxis growth system. J. Lond. Math. Soc. 74(2), 453–474 (2006)
  • (17) Strehl, R., Sokolov, A., Kuzmin, D., Horstmann, D., Turek, S.: A positivity-preserving finite element method for chemotaxis problems in 3D. J. Comput. Appl. Math. 239, 290–303 (2013)
  • (18) Strehl, R., Sokolov, A., Kuzmin, D., Turek, S.: ”A flux-corrected finite element method for chemotaxis problems. Comput. Meth. Appl. Math. 10(2), 219–232 (2010)
  • (19) Akhmouch, M., Benzakour Amine, M.: Semi-implicit finite volume schemes for a chemotaxis-growth model. Indag. Math. 27(3), 702–720 (2016)
  • (20) Horstmann, D.: From 1970 until now: the Keller-Segel model in chemotaxis and its consequences i. Jahresber. Dtsch. Math.-Ver. 105, 103–165 (2003)
  • (21) Hillen, T., Painter, K.J.: A user’s guide to PDE models for chemotaxis. J. Math. Biol. 58, 183–217 (2009)
  • (22) Perthame, B.: Transport Equations in Biology. Birkhäuser, Basel (2007)
  • (23) Winkler, M.: Boundedness in the higher-dimensional parabolic-parabolic chemotaxis system with logistic source. Commun. Partial Diff. Eqn. 35, 1516–1537 (2010)
  • (24) Hillen, T., Painter, K.: Global existence for a parabolic chemotaxis model with prevention of overcrowding. Adv. Appl. Math., 26(4), 280–301 (2001)
  • (25) Painter, K. J., Hillen, T.: Volume-filling and quorum-sensing in models for chemosensitive movement. Can. Appl. Math. Quart, 10(4), 501–543 (2002)
  • (26) Andreianov, B., Bendahmane, M., Saad, M.: Finite volume methods for degenerate chemotaxis model. J. Comput. Appl. Math. 235 (14), 4015–4031 (2011)
  • (27) Chamoun, G., Saad, M., Talhouk, R.: Monotone combined edge finite volume-finite element scheme for anisotropic Keller-Segel model. Numer. Methods Partial Differential Equation 30 (3), 1030–1065 (2014)
  • (28) Saito, N.: Conservative upwind finite-element method for a simplified Keller-Segel system modelling chemotaxis, IMA J. Numer. Anal. 27, 332–365 (2007)
  • (29) Zhou, G., Saito, N.: Finite volume methods for a Keller-Segel system: discrete energy, error estimates and numerical blow-up analysis. Numerische Mathematik, 135(1), 265–311 (2017)
  • (30) Ibrahim, M., Saad, M.: On the efficacy of a control volume finite element method for the capture of patterns for a volume-filling chemotaxis model. Comput. Math. Appl. 68(9), 1032–1051 (2014)
  • (31) Filbet, F.: A finite volume scheme for the Patlak-Keller-Segel chemotaxis model. Numer. Math. 104(4), 457–488 (2006)
  • (32) Lagaris, I. E., Likas, A., Fotiadis, D. I.: Artificial neural networks for solving ordinary and partial differential equations. IEEE Trans. Neural Netw. 9(5), 987–1000 (1998)
  • (33) Lagaris, I. E., Likas, A. C., Papageorgiou, D. G.: Neural-network methods for boundary value problems with irregular boundaries. IEEE Trans. Neural Netw. 11(5), 1041–1049 (2000)
  • (34)

    Shirvany, Y., Hayati, M., Moradian, R.: Multilayer perceptron neural networks with novel unsupervised training method for numerical solution of the partial differential equations. Appl. Soft Comput. 9(1), 20–29 (2009)

  • (35) Mai-Duy, N., Tran-Cong, T.: Numerical solution of differential equations using multiquadric radial basis function networks. Neural Netw. 14(2), 185–199 (2001)
  • (36) Aminataei, A., Mazarei, M. M.: Numerical solution of Poisson’s equation using radial basis function networks on the polar coordinate. Comput. Math. with Appl. 56(11), 2887–2895 (2008)
  • (37) Jianyu, L., Siwei, L., Yingjian, Q., Yaping, H.: Numerical solution of elliptic partial differential equation using radial basis function neural networks. Neural Netw. 16, 729–734 (2003)
  • (38) Sun, H., Hou, M., Yang, Y., Zhang, T., Weng, F., Han, F.: Solving partial differential equation based on Bernstein neural network and extreme learning machine algorithm. Neural Process. Lett. 50(2), 1153–1172 (2019)
  • (39) Mall, S., Chakraverty, S.: Single layer Chebyshev neural network model for solving elliptic partial differential equations. Neural Process. Lett. 45(3), 825–840 (2017)
  • (40) Rudd, K., Ferrari, S.: A constrained integration (CINT) approach to solving partial differential equations using artificial neural networks. Neurocomputing. 155, 277–285 (2015)
  • (41) Beidokhti, R. S., Malek, A.: Solving initial-boundary value problems for systems of partial differential equations using neural networks and optimization techniques. J. Franklin. Inst. 346(9), 898–913 (2009)
  • (42)

    Sirignano, J., Spiliopoulos, K.: DGM: A deep learning algorithm for solving partial differential equations. J. Comput. Phys. 375, 1339–1364 (2018)

  • (43) Weinan, E., Han, J., Jentzen, A.: Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations. Commun. Math. Stat. 5(4), 349–380 (2017)
  • (44) Eymard, R., Gallouët, T., Herbin, R.: Finite volume methods. In: P. G. Ciarlet and J. L. Lions(eds.), Handbook of numerical analysis volume VII, pp. 713–1020, North-Holland (2000)
  • (45) Spalding, D. B.: A novel finite difference formulation for differential expressions involving both first and second derivatives. Internat. J. Numer. Methods Engrg. 4(4), 551–559 (1972)
  • (46) Bessemoulin-Chatard, M., Chainais-Hillairet, C., Filbet, F.: On discrete functional inequalities for some finite volume schemes. IMA J. Numer. Anal. 35(3), 1125–1149 (2015)
  • (47) Chainais-Hillairet, C., Liu, J.-G., Peng, Y.-J: Finite volume scheme for multi-dimensional drift-diffusion equations and convergence analysis. Math. Mod. Numer. Anal. 37, 319–338 (2003)
  • (48) Yu, H., Wilamowski, B. M.: Levenberg-Marquardt training. Industrial electronics handbook, 5(12), pp. 1–12, (2011)
  • (49)

    Tan, H. H., Lim, K. H.: Review of second-order optimization techniques in artificial neural networks backpropagation. In: IOP Conference Series: Materials Science and Engineering (Vol. 495, No. 1, p. 012003). IOP Publishing (2019)

  • (50) Woodward, D. E., Tyson, R., Myerscough, M. R., Murray, J. D., Budrene, E. O., Berg, H. C.: Spatio-temporal patterns generated by Salmonella typhimurium. Biophysical journal, 68(5), 2181–2189 (1995)
  • (51) Thar, R., Kühl, M.: Complex pattern formation of marine gradient bacteria explained by a simple computer model. FEMS microbiology letters, 246(1), 75–79 (2005)
  • (52) Budrene, E.O., Berg, H.C.: Dynamics of formation of symmetrical patterns of chemotactic bacteria. Nature 376, 49–53 (1995)