Analytical Solutions to the Chavy-Waddy–Kolokolnikov Model of Bacterial Aggregates in Phototaxis by Three Integration Schemes

Abstract

In this work, we find analytical solutions to the Chavy-Waddy–Kolokolnikov equation, a continuum approximation for modeling aggregate formation in bacteria moving toward the light, also known as phototaxis. We used three methods to obtain the solutions, the generalized Kudryashov method, the $e^{-R\left(\xi \right)}$-expansion, and exponential function methods, all of them being very efficient for finding traveling wave-like solutions. Findings can be classified into the case where the nonlinear term can be considered a small perturbation of the linear case and the regime of instability and pattern formation. Standing waves and traveling fronts were also found among the physically interesting cases, in addition to recovering stationary spike-like solutions.

1. Introduction

Bacteria have evolved various mechanisms to adapt and survive in their environments, such as migrating towards regions with higher nutrient concentrations or better living conditions. Chemotaxis is a common adaptation where bacteria sense and respond to chemical gradients, moving towards a region with a higher concentration of a particular substance. Another adaptation is phototaxis, which involves the movement of photosynthetic, motile organisms towards light. Both chemotaxis and phototaxis play important roles in evolutionary and ecological processes. Chemotaxis has been extensively studied in biology and mathematics, with the Keller–Segel equation being one of the earliest and most well-known models [1]. This equation consists of a reaction–diffusion–advection-like equation for bacterial density containing a function for chemotactic sensitivity, another function for the production and death of individuals, and a cross-diffusion term that couples with the concentration of the chemical signal that has its kinetics [2]. On the other hand, from a biological perspective, it has been verified that for the successful realization of phototaxia, the presence of both photoreceptors and pili is crucial, as they play a key role in facilitating its progression [3], which proved to be fundamental in agent model simulations [4]. From mathematical modeling, there is a series of papers, where D. Levy et al. [5,6,7,8,9,10] proposed some models to describe how phototaxis bacteria behave based on some basic features extracted from observations and experiments. The tools range from stochastic equations, particle models, kinetic models, to master equations in terms of probabilities and partial differential equations (PDEs). Group dynamics is a crucial feature in this system and was encoded through an internal degree of freedom called excitation [6,9]. Some models show additional internal variables involved with excitation, such as rotation; it was found that the sensitivity to perform phototaxis decreases if there is no rotation of the colony [11]. Among Levy’s papers, there are two relevant PDE models, the first being a reaction–diffusion equation system for bacterial concentration, excitation, and substrate memory derived as a continuous limit of a stochastic model [5,6]. This resembles the chemotactic model, which has also been adapted to analyze phototaxis [12,13]. The second model involves local interactions by means of a proposed system of master equations for the probability of finding bacteria in a particular state. Such a system includes reaction–diffusion, persistence, and sticking terms [9,10]. From the latter, Chavy-Waddy and Kolokolnikov (CWK) proposed a simplified system of equations for the probability that the bacteria move and obtained a fourth-order nonlinear partial differential equation, only depending on one parameter that combines the probabilities of moving to either side, staying, or changing direction according to the sensing distance [14]. The resulting model is of swarming type for bacterial aggregation, like the Cahn–Hillard equation [15]. The stationary solution coincides with a state of particle aggregation, for which Taranets and Chugunova [16] studied the rate of convergence and the existence of non-negative solutions. Recently, the physical characteristics of the bacteria, such as their shape and the way the flaps work, are also being explored for increasingly accurate models of phototaxis [17]. The CWK equation includes a reverse diffusion term, a fourth-order term related to a long-range effect term, and a nonlinear term with a unique parameter that considers the aggregate extent and whose value determines the instability region for structure formation.

In this paper, we solve the CWK equation for propagating non-deformable pulses, employing three generalizations of Kudryashov’s method. The Kudryashov method is highly efficient for finding exact solutions to nonlinear differential equations. It has a wide range of applications, from physics, engineering, mathematics, and biology. Particularly in biology, it was used to delve into nonlinear phenomena, such as the study of HIV-1 infection [18], population dynamics [19], etc. These methods are suitable for describing soliton-like traveling wave solutions that have a clear biophysical interpretation in the present model.

The article’s structure is the following. We present the CWK phototaxis model that will be solved in Section 2. Section 3 provides an overview of the methods employed, including the solutions found in each case and some graphical representations of important cases. In Section 4, we highlight the importance of our solutions and discuss the overall results. We also present five appendices with the largest expressions and additional results. Our findings provide a collection of exact solutions for studying phototaxis from the reduced CWK one-dimensional model.

2. Mathematical 1D Model of Aggregation in Bacterial Colonies

As mentioned previously, Chavy-Waddy and Kolokolnikov proposed in [14] a nonlinear parabolic fourth-order partial differential equation for modeling the movement of phototaxic bacterial aggregates using the random models proposed by Levy et al. The CWK formula is as follows:

(1)$u_t=-u_{xx}-u_{xxxx}+\alpha {\left(\frac{u_x{u}_{xx}}{u}\right)}_x.$

The first two right-hand terms are similar to reverse diffusion and a long-range term. Reverse diffusion occurs when transport is towards zones where the concentration gradient is high, opposite to what happens in diffusion. This is the case, for instance, in the Cahn–Hilliard equation for the phase separation process [15,20]. The fourth-order term is sometimes associated with long-range terms where the influence of distant neighbors on the concentration at a given point is considered [21,22]. Both terms balance each other; while the inverse diffusion destabilizes the system, the fourth-order term stabilizes the higher Fourier modes. The third term is the one containing the nonlinearity and the only parameter of the model, $\alpha$, which controls the size of the aggregate and is given by

(2)$\alpha =\frac{c(2d+1){(d+1)}^2}{\left(c[1+d(d^2+2d+3)]-2a\right)},$

where the constants are given in terms of the simplified Levy’s master equation [9]. a is the jump rate at which the bacterium moves, preserving its orientation, c is the rate at which it moves after switching to a new orientation, and d is the bacterium’s sensing radius. The model given by Equation (1) is similar to the swarming models of biological aggregation based on attraction–repulsion forces [23,24]. However, it is a way simpler since it only involves a single equation with stationary finger-like solutions, as found in the experiments. The stationary finger-like solution found in [14] is as follows:

(3)$u(x,t)=A{\left[\sec \mathrm{h}\left(\frac{\sqrt{\alpha -1}}{2}x\right)\right]}^{\frac{2}{\alpha -1}},$

where A is a real normalization constant, and only depends on the value of $\alpha$. We note that $\alpha =1$ is a particular value; indeed, in [14], it is shown that $\alpha >1$ is necessary to obtain the steady state since it is obtained in the unstable regime when $c>2a/d$, where patterns can occur; see Appendix A. Some extreme cases satisfy this condition when the motion rate after changing orientation becomes very large $c\to \infty$. If the bacterium stops without changing orientation $a=0$, one obtains $\alpha \ge \left[{(1+d)}^2(1+2d)\right]/[1+d(3+2d+d^2)]$. From this, three cases follow depending on the value of the sensitivity distance. If $d=1$, then $\alpha \ge 12/7$; if $d\to \infty$, then $\alpha \ge 2$; if $d=0$, then $\alpha \ge 1$. In all these cases, the unstable threshold is fulfilled in general when $\alpha >1$.

The simplest case is when $\alpha =0$, where Equation (1) reduces to $u_t=-u_{xx}-u_{xxxx}$, which only has a contribution to the flux due to inverse diffusion and long-range terms. It is a linear equation that conventional methods can solve; nonetheless, the methods used in this paper also give additional solutions, so we will present them for the sake of completeness. Appendix B presents the solutions of the case $\alpha =0$ by the method of separation of variables and comparison with some of the solutions obtained here. The nonlinear equation is not straightforward to solve by the usual methods, so for non-zero $\alpha$ cases, solutions are obtained by the methods explored herein and include two regimes, $0<\alpha <1$ and $\alpha >1$, being the most relevant cases.

Finding the stationary solution of Equation (1) is accomplished in [14] by reducing the order of the equation and then studying the orbits of the system, which involves the following transformation $u(x,t)=e^{v(x,t)}$, and Equation (1) becomes

(4)$v_t=-v_x^2-v_{xx}-v_{xxxx}+(4\alpha -6)v_x^2{v}_{xx}+(\alpha -4)v_x{v}_{xxx}+(\alpha -3)v_{xx}^2+(\alpha -1)v_x^4.$

In this equation, there seem to be four particular values of $\alpha$. However, when changing to $z=v_x$, the corresponding equation has only two characteristic values for $\alpha =1,3$. Indeed, for $\alpha =3$, the equation for z can be easily integrated, as found in [16], where they also analyze the stability of some stationary state families of the CWK equation. Moreover, time-dependent solutions were presented, giving their convergence rate to the stationary case.

In the next section, we will present some stationary and time-dependent solution families, each obtained with the so-called Kudryashov method. Since Equation (1) is of the parabolic type, and it is well known that parabolic equations admit traveling wave-like solutions as in [25,26], we expect that the CWK equation also admits soliton-like solutions suitable to be found with the Kudryashov method. To do so, we use Equation (4) and make the following variable change $\xi =x-\omega t$, under the assumption of traveling-wave-like solutions, resulting in the following:

(5)$\omega v_{\xi }-v_{\xi }^2-v_{\xi \xi }-v_{\xi \xi \xi \xi }+(4\alpha -6)v_{\xi }^2{v}_{\xi \xi }+(\alpha -4)v_{\xi }{v}_{\xi \xi \xi }+(\alpha -3)v_{\xi \xi }^2+(\alpha -1)v_{\xi }^4=0,$

where $\omega$ is the constant wave velocity. Lastly, we substitute $\varphi \left(\xi \right)=v_{\xi }$ to obtain a simplified version of Equation (5):

(6)$\omega \varphi -{\varphi }^2-{\varphi }_{\xi }-{\varphi }_{\xi \xi \xi }+(4\alpha -6){\varphi }^2{\varphi }_{\xi }+(\alpha -4)\varphi {\varphi }_{\xi \xi }+(\alpha -3){\left({\varphi }_{\xi }\right)}^2+(\alpha -1){\varphi }^4=0.$

We will present in the next section the methods used to solve Equation (6), and hence Equation (1), and the collection of families that each technique will produce.

3. Methods and Solutions

In this section, we will describe briefly each of the proposed methods, their application, and the families of solutions obtained by means of them.

3.1. Brief Description of the Generalized Kudryashov Method

The purpose of this section is to present the algorithm of the generalized Kudryashov method for finding exact solutions of nonlinear evolution equations, such as Equation (6), consisting of the following steps:

• Step 1: 

  We assume the exact solutions to Equation (6) can be formulated as follows:

  (7)$\varphi \left(\xi \right)=\frac{{\sum }_{i=0}^N{a}_i{Q}^i\left(\xi \right)}{{\sum }_{j=0}^M{b}_i{Q}^j\left(\xi \right)},$

  where $a_i$ and $b_j$ are arbitrary constants with $a_N\ne 0,b_M\ne 0,$ and the function $Q\left(\xi \right)$ satisfies the next differential equation [27]:

  (8)$Q_{\xi }=Q^2-Q.$

  The relation between integers N and M can be established by considering the homogeneous balance between the higher-order derivatives and the nonlinear factors in Equation (6). In our case, $N=2$ and $M=1$.

• Step 2: 

  Next we substitute both $\varphi$, given in Equation (7), and its derivatives ${\varphi }_{\xi },{\varphi }_{\xi \xi },\dots ,$ in Equation (6) to obtain the polynomial equation:

  (9)$P\left(Q\left(\xi \right)\right)=0.$

• Step 3: 

  We select all the terms having the same algebraic power in Q from the polynomial Equation (9), setting them equal to zero, and obtain a system of algebraic equations with the following set of unknowns, $\{a_0,a_1,a_2,b_0,b_1,\omega \}$ depending on the value of $\alpha$. We use algebraic manipulation software such as Mathematica to solve the system with the model constraints, considering that $a_2\ne 0$ and $b_1\ne 0$ are also required.

• Step 4: 

  Using the previous results and considering Equation (7) together with Equation (8), we obtain the possible exact solutions of Equation (6) and consequently those of Equation (1).

  Due to the fact that the generalized Kudryashov method is defined by the rational form of finite series given by Equation (7), it provides a greater number of exact and more general solutions in an identical manner as the classical Kudryashov method, which is a significant advantage [28].

Solutions Obtained by the Generalized Kudryashov Technique

The system of nonlinear algebraic equations resulting from this method is shown in Appendix C in Equation (A27). Next, we present the solutions obtained for different values of the parameters. The first set of solutions is for $\alpha =0$.

Solution 1.

If $\alpha =0$ , we have $a_0=a_0$ , $a_1=-a_2-b_0$ , $a_2=a_2$ , $b_0=a_0$ , $b_1=-a_2$ , and $\omega =2$ , from which we obtain the solution

(10) $u_1(x,t)=cosh(2t-x)-sinh(2t-x)+1.$

Solution 2.

If $\alpha =0$ , we have $a_0=0$ , $a_1=-b_0$ , $a_2=a_2$ , $b_0=b_0$ , $b_1=-a_2$ , and $\omega =-2$ , from which we obtain the solution

(11) $u_2(x,t)=-sinh(2t+x)+cosh(2t+x)+1.$

Solution 3.

If $\alpha =0$ , we have $a_0=0$ , $a_1=0$ , $a_2=2b_0$ , $b_0=b_0$ , $b_1=-a_2$ , and $\omega =-10$ , from which we obtain the solution

(12) $u_3(x,t)=-sinh(20t+2x)+cosh(20t+2x)-1.$

Solution 4.

If $\alpha =0$ , we have $a_0=2b_0$ , $a_1=-4b_0$ , $a_2=2b_0$ , $b_0=b_0$ , $b_1=-a_2$ , and $\omega =10$ , from which we obtain the solution

(13) $u_4(x,t)=sinh(20t-2x)-cosh(20t-2x)+1.$

In all cases, non-stationary waves propagating in different directions were obtained. To illustrate, let us consider Solution 4. In Figure 1, we show the plot of traveling wave solution $u_4(x,t)$.

Next we present the solutions for $\alpha \ne 0$.

Solution 5.

If $\alpha =1$ , we have $a_0=\frac{b_0(a_1{b}_1-a_2{b}_0)}{b_1^2}$ , $a_1=a_1$ , $a_2=a_2$ , $b_0=b_0$ , $b_1=b_1$ , and $\omega =0$ , from which we obtain the solution

(14) $u_5\left(x\right)=cosh\left(\frac{1}{4}x(\beta +x)\right)-sinh\left(\frac{1}{4}x(\beta +x)\right),\beta =\frac{4a_2{b}_0-4b_1(2a_2+a_1)}{b_1^2}.$

Solution 6.

If $\alpha >1$ , we have $a_0=\frac{1}{4}(-2a_1-a_2)$ , $a_1=a_1$ , $a_2=a_2$ , $b_0=-a_1-\frac{1}{2}{a}_2$ , $b_1=-a_2$ , and $\omega =0$ , from which we obtain the solutions

(15) $u_6\left(x\right)=Ccosh\left(\gamma x\right),C\ne 0and\gamma =\pm \frac{1}{\sqrt{\alpha -1}}.$

Solution 7.

If $\alpha >1$ , we have $a_0=\frac{1}{2}{a}_2$ , $a_1=-a_2$ , $a_2=a_2$ , $b_0=-a_1-\frac{1}{2}{a}_2$ , $b_1=-a_2$ , and $\omega =0$ , from which we obtain the solutions

(16) $u_7\left(x\right)=Csinh\left(\gamma x\right),C\ne 0and\gamma =\pm \frac{1}{\sqrt{\alpha -1}}.$

Solution 8.

If $\alpha >1$ , we have $a_0=\frac{1}{4}(-2a_1-a_2)$ , $a_1=a_1$ , $a_2=a_2$ , $b_0=\frac{1}{4}(\alpha -1)(2a_1+a_2)$ , $b_1=\frac{1}{2}(\alpha -1)a_2$ , and $\omega =0$ , from which we obtain the solutions

(17) $u_8\left(x\right)=C{\mathit{sec}\mathrm{h}}^m\left(\gamma x\right),C\ne 0,m=\frac{2}{\alpha -1}and\gamma =\frac{\sqrt{\alpha -1}}{2}.$

Here only stationary solutions with $\omega =0$ were obtained, recovering especially the finger-like distribution from [14] in Solution 8. Figure 2 shows that the bell-shaped curve’s distribution becomes progressively wider as the value of alpha increases. This suggests that the bacteria exhibit a preference for being farther away. In other words, as alpha increases, the bacteria tend to distribute themselves over a larger area, indicating a lower degree of clustering.

Earlier, we mentioned that when $\alpha =3$, Equation (4) for $v_x$ can be directly integrated. For this special case, the following solutions were found.

Solution 9.

If $\alpha =3$ , we have $a_0=\frac{a_2}{12-6\sqrt{3}}$ , $a_1=-\frac{1}{3}(\sqrt{3}+3)a_2$ , $a_2=a_2$ , $b_0=-a_1-\frac{1}{2}{a}_2$ , $b_1=-a_2$ , and $\omega =\frac{1}{3\sqrt{3}}$ , from which we obtain the solution

(18) $u_9(x,t)=-\frac{e^{\frac{x}{\sqrt{3}}-\frac{t}{9}}\left(e^{\frac{t}{3\sqrt{3}}}+e^x\right)}{\left(2\sqrt{3}-3\right)e^{\frac{t}{3\sqrt{3}}}+\left(2\sqrt{3}+3\right)e^x}$

Solution 10.

If $\alpha =3$ , we have $a_0=\frac{a_2}{6(\sqrt{3}+2)}$ , $a_1=\frac{1}{3}(\sqrt{3}-3)a_2$ , $a_2=a_2$ , $b_0=-a_1-\frac{1}{2}{a}_2$ , $b_1=-a_2$ , $\beta \in \mathbb{R}$ , $h=\pm 1$ , and $\omega =-\frac{1}{3\sqrt{3}}$ , from which we obtain the solution

(19) $\begin{array}{cc}\hfill u_{10}(x,t)& =\beta \left[sinh\left(\frac{t}{3\sqrt{3}}+hx\right)+cosh\left(\frac{t}{3\sqrt{3}}+hx\right)+1\right]\left[\sqrt{3}sinh\left(\frac{t}{3\sqrt{3}}+hx\right)\right.\hfill \\ \hfill & {\left.+\sqrt{3}cosh\left(\frac{t}{3\sqrt{3}}+hx\right)+7\sqrt{3}+12\right]}^{\frac{1}{2\sqrt{3}}-\frac{1}{2}}\left[\left(2\sqrt{3}+3\right)\left(sinh\left(\frac{t}{3\sqrt{3}}+hx\right)\right.\right.\hfill \\ \hfill & {\left.\left.+cosh\left(\frac{t}{3\sqrt{3}}+hx\right)\right)+26\sqrt{3}+45\right]}^{-\frac{1}{2}-\frac{1}{2\sqrt{3}}}\left[cosh\left(\frac{t}{9}+\frac{hx}{\sqrt{3}}\right)-sinh\left(\frac{t}{9}+\frac{hx}{\sqrt{3}}\right)\right]\hfill \end{array}$

Solution 11.

If $\alpha =3$ , we have $a_0=\pm \frac{1}{6}(\sqrt{3}\pm 2)a_2$ , $a_1=\mp \frac{1}{3}(\sqrt{3}\pm 3)a_2$ , $a_2=a_2$ , $b_0$ = $\frac{1}{4}(\alpha -1)(2a_1+a_2)$ , $b_1=\frac{1}{2}(\alpha -1)a_2$ , and $\omega =\mp \frac{1}{3\sqrt{3}}$ , from which we obtain the solution

(20) $u_{{11}_{\mp }}(x,t)=\mp \left(3tanh\left(\frac{1}{18}\left(\sqrt{3}t\pm 9x\right)\right)+2\sqrt{3}\right)e^{\mp \frac{x}{\sqrt{3}}-\frac{t}{9}}.$

The above solutions exhibit similar behaviors at different scales. Although they have exponential growth, near zero, there is a small propagating pulse. To illustrate this, we consider the solution $u_{{11}_{-}}(x,t)$. In Figure 3, we present the graph of the traveling pulse of Solution 11; notice how it moves to the right.

Since the method allows specific solutions for particular values of the parameters, here we present the stationary solution for $\alpha =5$, being a particular case for $\alpha >1$.

Solution 12.

If $\alpha =5$ , we have $a_0=\frac{1}{2}{a}_2$ , $a_1=-a_2$ , $a_2=a_2$ , $b_0=\frac{1}{4}(\alpha -1)(2a_1+a_2)$ , $b_1=\frac{1}{2}(\alpha -1)a_2$ , and $\omega =0$ , from which we obtain the solution

(21) $u_{12}\left(x\right)=\frac{sinh\left(\frac{x}{2}\right)+cosh\left(\frac{x}{2}\right)}{\sqrt{-2{sinh}^2\left(x\right)-2sinh\left(x\right)cosh\left(x\right)}}.$

3.2. Brief Description of the $e^{-R\left(\xi \right)}$-Expansion Method

Among the methods for finding analytical solutions to nonlinear equations is the so-called $e^{-R\left(\xi \right)}$-expansion, which has been used to find solitary wave-like solutions in some fluid problems [29,30].

• Step 1: 

  The $e^{-R\left(\xi \right)}$-expansion method assumes that the solution of Equation (6) is expressed as

  (22)$\varphi \left(\xi \right)=\sum _{i=0}^N{a}_i{\left(e^{-R\left(\xi \right)}\right)}^i$

  where $a_i$ is an arbitrary constant with $a_N\ne 0$, and the function R satisfies the following differential equation [31]:

  (23)$R_{\xi }=\lambda +\mu e^{R\left(\xi \right)}+e^{-R\left(\xi \right)}.$

  Consequently, the function R can be given by

  (24)$R\left(\xi \right)=\left\{\begin{array}{cc}ln\left(-\frac{\sqrt{{\lambda }^2-4\mu }tanh\left(\frac{1}{2}(\xi +A)\sqrt{{\lambda }^2-4\mu }\right)}{2\mu }-\frac{\lambda }{2\mu }\right)\hfill & \mathrm{if}{\lambda }^2-4\mu >0,\mu \ne 0\hfill \\ ln\left(-\frac{\sqrt{{\lambda }^2-4\mu }coth\left(\frac{1}{2}(\xi +A)\sqrt{{\lambda }^2-4\mu }\right)}{2\mu }-\frac{\lambda }{2\mu }\right)\hfill & \mathrm{if}{\lambda }^2-4\mu >0,\mu \ne 0\hfill \\ ln\left(\frac{\sqrt{4\mu -{\lambda }^2}tan\left(\frac{1}{2}(\xi +A)\sqrt{4\mu -{\lambda }^2}\right)}{2\mu }-\frac{\lambda }{2\mu }\right)\hfill & \mathrm{if}{\lambda }^2-4\mu <0,\mu \ne 0\hfill \\ ln\left(\frac{\sqrt{4\mu -{\lambda }^2}cot\left(\frac{1}{2}(\xi +A)\sqrt{4\mu -{\lambda }^2}\right)}{2\mu }-\frac{\lambda }{2\mu }\right)\hfill & \mathrm{if}{\lambda }^2-4\mu <0,\mu \ne 0\hfill \\ -ln\left(\frac{\lambda }{e^{\lambda (\xi +A)}-1}\right)\hfill & \mathrm{if}\mu =0,\lambda >0\hfill \\ ln\left(-\frac{2\left(\lambda \right(\xi +A)+2)}{{\lambda }^2(\xi +A)}\right)\hfill & \mathrm{if}\lambda \ne 0,{\lambda }^2-4\mu =0\hfill \\ ln(\xi +A)\hfill & \mathrm{if}\mu =0,\lambda =0.\hfill \end{array}\right.$

  As previously said, to compute the positive integer N, consider the homogeneous balance between the higher-order derivatives and the nonlinear parts in Equation (6). In this case, $N=1$.

• Step 2: 

  In this method we consider $\varphi$ given in Equation (22) and the necessary derivatives ${\varphi }_{\xi },{\varphi }_{\xi \xi },\dots ,$ then we substitute them into Equation (6) to obtain the following polynomial equation:

  (25)$P\left(e^{-R\left(\xi \right)}\right)=0.$

• Step 3: 

  We select from the polynomial Equation (25) all terms having the same algebraic power of $e^{-R\left(\xi \right)}$, set them equal to zero, and obtain a system of algebraic equations with the set of unknowns $\{a_0,a_1,\omega \}$ depending on $\alpha$. In the same way as the previous method, we use Mathematica to solve the system with its natural constraints, assuming $a_1\ne 0$.

• Step 4: 

  With the results obtained in the previous step and taking Equation (22) along with Equation (23), we obtain the possible exact solutions of Equation (6), and hence those of Equation (1).

Solutions Found by the $e^{-R\left(\xi \right)}$-Expansion Method

The resulting nonlinear algebraic system of equations resulting from this method is presented in Equation (A28) of Appendix D; here we present the solutions obtained by this method. First for the case $\alpha =0$:

Solution 13.

If $\alpha =0$ , we have $a_0=0$ , $a_1=1$ , $\lambda >0$ , $\mu =0$ , and $\omega =-{\lambda }^3-\lambda$ , from which we obtain the solution

(26) $u_{13}(x,t)=e^{-\lambda \left(\left({\lambda }^3+\lambda \right)t+x\right)}-e^{A\lambda },A\in \mathbb{R}.$

Solution 14.

If $\alpha =0$ , we have $a_0>0$ , $a_1=1$ , $\lambda =a_0$ , $\mu =a_0(\lambda -a_0)$ , and $\omega =6a_0{\lambda }^2-12a_0^2\lambda +8a_0^3+2a_0-{\lambda }^3-\lambda$ , from which we obtain the solution

(27) $u_{14}(x,t)=-sinh\left(a_0(x+A)-a_0^{2}t(a_0^2+1)\right)-cosh\left(a_0(x+A)-a_0^{2}t(a_0^2+1)\right)+1,A\in \mathbb{R}.$

The first two solutions are similar to $u_4(x,t)$, with wavefronts propagating to one side. Figure 4 shows, with fixed values of the parameters, a constant unit density over time.

Solution 15.

If $\alpha =0$ , we have $a_0=a_0$ , $a_1=2$ , $\lambda =a_0$ , $\mu =\frac{1}{4}\left(a_0^2+1\right)$ , and $\omega =0$ , from which we obtain the solution

(28) $u_{15}\left(x\right)={\left(sin\left(\frac{x+A}{2}\right)-a_{0}cos\left(\frac{x+A}{2}\right)\right)}^2,A\in \mathbb{R}.$

The solution $u_{15}\left(x\right)$ is an oscillatory function; evidently, the constants A and $a_0$ are the phase and amplitude of a standing wave. Notably, this solution arises only from the reverse diffusion and the fourth-order term. Thus, although each maximum represents the bacterial concentration, this distribution is preserved from the beginning of the process.

This method made it possible to find two stationary solutions for $\alpha <1$. While this clearly does not correspond to the region of instability, and therefore we cannot expect the formation of aggregates, we can think of $\alpha$ as a perturbation parameter in an intermediate region between reverse diffusion alone and pattern formation.

Solution 16.

If $\alpha =1-\frac{2}{a_1}$, we have $a_0=a_0$, $a_1>0$, $\lambda =\frac{2a_0}{a_1}$, $\mu =\frac{2a_0^2+a_1}{2a_1^2}$, and $\omega =0$, from which we obtain the solution

(29)$u_{16}\left(x\right)={\left(\sqrt{2}sin\left(\frac{x+A}{\sqrt{2a_1}}\right)-2a_0\sqrt{\frac{1}{a_1}}cos\left(\frac{x+A}{\sqrt{2a_1}}\right)\right)}^{a_1},A\in \mathbb{R}.$

Particularly, it makes sense for $a_1\ge 2$, which implies $\alpha <1$.

Solution 17.

If $\alpha =\frac{1}{3}$ , we have $a_0=a_0$ , $a_1=3$ , $\lambda =\frac{2a_0}{3}$ , $\mu =\frac{1}{18}\left(2a_0^2+3\right)$ , and $\omega =0$ , from which we obtain the solution

(30) $u_{17}\left(x\right)={\left(\sqrt{6}sin\left(\frac{x+A}{\sqrt{6}}\right)-2a_{0}cos\left(\frac{x+A}{\sqrt{6}}\right)\right)}^3,A\in \mathbb{R}.$

Interestingly, the solutions in this regime are also oscillatory, with Solution 16 having the same structure as $u_{17}$ for even values of $a_1$. However, for odd $a_1$, negative values occur, and the parity of $a_1$ must be considered to interpret $u_{16}$ as a distribution.

Figure 5 illustrates the solution with different values of $a_1$, which represents stationary distributions of bacteria aggregates. As the value of $a_1$ increases, the number of bacterial aggregates decreases while the amplitude of each curve increases. This leads to a higher density of bacteria within each curve.

Finally, in the region of structure formation, when $\alpha >1$, we obtain five stationary solutions, which also depend on method parameters.

Solution 18.

If $\alpha >1$ , we have $a_0=a_0$ , $a_1=-\frac{2}{\alpha -1}$ , $\lambda =\frac{2a_0}{a_1}$ , $\mu =\frac{2a_0^2+a_1}{2a_1^2}$ , and $\omega =0$ , from which we obtain the solution

(31) $u_{18}\left(x\right)={\left[\frac{e^{-\sqrt{\alpha -1}}-e^{\sqrt{\alpha -1}}}{\sqrt{\alpha -1}}+a_0{e}^{-\frac{1}{2}\sqrt{\alpha -1}(A+x)}\left(1+e^{\sqrt{\alpha -1}(A+x)}\right)\right]}^{-\frac{2}{\alpha -1}},A\in \mathbb{R}.$

Solution 19.

If $\alpha >1$ , we have $a_0=a_0$ , $a_1=1$ , $\lambda =2a_0$ , $\mu =a_0^2-\frac{1}{1-\alpha }$ , and $\omega =0$ , from which we obtain the solution

(32) $u_{19}\left(x\right)=\frac{1}{\sqrt{\alpha -1}}sinh\left(\frac{x+A}{\sqrt{\alpha -1}}\right)\mp a_{0}cosh\left(\frac{x+A}{\sqrt{\alpha -1}}\right),A\in \mathbb{R}.$

Solution 20.

If $\alpha =1+\frac{1}{a_0^2}>1$ , we have $a_0>0$ , $a_1=1$ , $\lambda =2a_0$ , $\mu =a_0^2-\frac{1}{1-\alpha }$ , and $\omega =0$ , from which we obtain the solution

(33) $u_{20}\left(x\right)=e^{-a_{0}x}\left(1-e^{2a_0(A+x)}\right),A\in \mathbb{R}.$

Solution 21.

If $\alpha =1-\frac{2}{a_1}$, we have $a_0<0$, $a_1=-2a_0^2$, $\lambda =\frac{2a_0}{a_1}$, $\mu =\frac{2a_0^2+a_1}{2a_1^2}$, and $\omega =0$, from which we obtain the solution

(34)$u_{21}\left(x\right)=e^{-a_{0}x}{\left(1-e^{-\frac{x+A}{a_0}}\right)}^{2a_0},A\in \mathbb{R}.$

Note that given the restriction for $a_1$, $\alpha =1+\frac{1}{a_1^2}>1$.

Solution 22.

If $\alpha =1-\frac{2}{a_1}$, we have $a_0=a_0$, $a_1<0$, $\lambda =\frac{2a_0}{a_1}$, $\mu =\frac{2a_0^2+a_1}{2a_1^2}$, and $\omega =0$, from which we obtain the solution

(35)$u_{22}\left(x\right)={\left(\sqrt{-\frac{2}{a_1}}{a}_{1}cosh\left(\frac{x+A}{\sqrt{-2a_1}}\right)+2a_{0}sinh\left(\frac{x+A}{\sqrt{-2a_1}}\right)\right)}^{a_1},A\in \mathbb{R}.$

Note that $\alpha >1$, given that $a_1<0$.

Some of the solutions are exponential, but there are also spike solutions, as in Solution 18.

Figure 6 illustrates, for $\alpha =5$, that the distribution of the bacterial population takes a bell-shaped curve in the steady state. In this scenario, the bacteria form an aggregate, and the density of the aggregate is determined by the value of $a_0$. When this value increases, the distance required for the bacteria to join and form an aggregate also increases. As a result, the bacteria are farther apart from each other, leading to a lower overall density of aggregates.

3.3. Brief Description of the Modified Exponential Function Method

The Exp-method was introduced to find solitary, compact, and periodic solutions of nonlinear wave-like equations [32]. It has been applied, for instance, to obtain soliton-type solutions for the Allen–Cahn equation, a reaction–diffusion equation describing phase separation in multi-alloy systems, and plasma dynamics [33]. The algorithm of the method is given below.

• Step 1: 

  We assume the exact solutions to Equation (6) can also be formulated as follows:

  (36)$\varphi \left(\xi \right)=\frac{{\sum }_{i=0}^N{a}_i{\left(e^{-Q\left(\xi \right)}\right)}^i}{{\sum }_{j=0}^M{b}_i{\left(e^{-Q\left(\xi \right)}\right)}^j},$

  where the $a_i$ and $b_j$ are arbitrary constants with $a_N\ne 0,b_M\ne 0,$ and the function Q satisfies the differential equation [32,33]:

  (37)$Q_{\xi }=\lambda +\mu e^{Q\left(\xi \right)}+e^{-Q\left(\xi \right)}.$

  Consequently, the function Q satisfies the same differential equation given in Equation (24).

  The integers N and M that appear in this method can be determined in the same way as before by considering the homogeneous balance between the higher-order derivatives and the nonlinear factors in Equation (6). In this case, $N=2$ and $M=1$.

  The second, third, and fourth steps of the current procedure are identical to those outlined in Section 3.2.

Solutions Found by the Modified Exponential Function Method

The nonlinear algebraic system of equations necessary to obtain solutions according to the exponential function method can be seen in Equation (A29) of Appendix E. Next we show the solutions we obtain by this method. For $\alpha =0$, we find three stationary and three traveling solutions below.

Solution 23.

If $\alpha =0$, $a_0=0$, $a_1=a_1$, $a_2=2b_1$, $b_0=0$, $b_1\ne 0$, $\lambda =\frac{a_1}{b_1}$, $\mu =\frac{a_1^2+b_1^2}{4b_1^2}$ and $\omega =0$, from which we obtain the solution

(38) $u_{23}\left(x\right)={\left(a_{1}cos\left(\frac{x+A}{2}\right)-b_{1}sin\left(\frac{x+A}{2}\right)\right)}^2,A\in \mathbb{R}.$

Solution 24.

If $\alpha =0$, $a_0\ne 0$, $a_1=0$, $a_2=-\frac{a_0}{\mu }$, $b_0=0$, $b_1=-\frac{a_0}{\mu }$, $\lambda =\mp \sqrt{4\mu -1}$, $\mu \ne 0$ and $\omega =0$, from which we obtain the solution

(39) $u_{24}\left(x\right)=\frac{1}{2}\left(\sqrt{4\mu -1}cos(x+A)\pm sin(x+A)+\sqrt{4\mu -1}\right),A\in \mathbb{R}.$

Solution 25.

If $\alpha =0$, $a_0\ne 0$, $a_1=0$, $a_2=-\frac{a_0}{\mu }$, $b_0=0$, $b_1=-\frac{a_0}{\mu }$, $\lambda =\sqrt{4\mu -1}$, $\mu \ne 0$ and $\omega =0$, from which we obtain the solution

(40) $u_{25}\left(x\right)=\frac{1}{2}\left(\sqrt{4\mu -1}cos(x+A)-sin(x+A)+\sqrt{4\mu -1}\right),A\in \mathbb{R}.$

Solution 26.

If $\alpha =0$, $a_0=0$, $a_1=a_1$, $a_2=b_1$, $b_0=a_1$, $b_1\ne 0$, $\lambda \ne 0$, $\mu =0$ and $\omega =-{\lambda }^3-\lambda$, from which we obtain the solution

(41) $u_{26}(x,t)=sinh\left(A\lambda \right)+cosh\left(A\lambda \right)+sinh\left({\lambda }^{4}t+{\lambda }^{2}t+\lambda x\right)-cosh\left({\lambda }^{4}t+{\lambda }^{2}t+\lambda x\right),A\in \mathbb{R}.$

Solution 27.

If $\alpha =0$, $a_0=0$, $a_1=0$, $a_2=b_1$, $b_0=\frac{b_1\lambda }{2}$, $b_1\ne 0$, $\lambda \ne 0$, $\mu =0$ and $\omega =-2\left(4{\lambda }^3+\lambda \right)$, from which we obtain the solution

(42) $u_{27}(x,t)=sinh\left(2A\lambda \right)+cosh\left(2A\lambda \right)+sinh\left(16{\lambda }^{4}t+4{\lambda }^{2}t+2\lambda x\right)-cosh\left(16{\lambda }^{4}t+4{\lambda }^{2}t+2\lambda x\right),A\in \mathbb{R}.$

Solution 28.

If $\alpha =0$, $a_0=0$, $a_1\ne 0$, $a_2=b_1$, $b_0=0$, $b_1\ne 0$, $\lambda =\frac{a_1}{b_1}$, $\mu =0$ and $\omega =\frac{a_1\left(a_1^2+b_1^2\right)}{b_1^3}$, from which we obtain the solution

(43) $u_{28}(x,t)=sinh\left(-\lambda A+{\lambda }^2({\lambda }^2+1)t-\lambda x\right)-cosh\left(-\lambda A+{\lambda }^2({\lambda }^2+1)t-\lambda x\right)+1,A\in \mathbb{R}.$

In this case, we have three standing wave-like solutions and three traveling wavefront-like solutions, close to those obtained with previous methods. Additionally, these solutions also depend on the parameters of the method.

Seven stationary solutions were found for the pre-pattern formation region, $\alpha <1$, all oscillatory functions, which are presented next.

Solution 29.

If $\alpha <1$, $a_0=-\frac{b_1}{8}$, $a_1=0$, $a_2=b_1$, $b_0=0$, $b_1\ne 0$, $\lambda =0$, $\mu =\frac{1}{8}$ and $\omega =0$, from which we obtain the solution

(44) $u_{29}\left(x\right)=Csin\left(\frac{x+B}{\sqrt{1-\alpha }}\right),C\ne 0,B\in \mathbb{R}.$

Solution 30.

If $\alpha <1$, $a_0=0$, $a_1=a_1$, $a_2=b_1$, $b_0=0$, $b_1\ne 0$, $\lambda =\frac{2a_1}{b_1}$, $\mu =\frac{\alpha a_1^2-a_1^2-b_1^2}{(\alpha -1)b_1^2}$ and $\omega =0$, from which we obtain the solution

(45) $u_{30}\left(x\right)=a_{1}cos\left(\sqrt{\frac{1}{1-\alpha }}(x+A)\right)-\sqrt{\frac{1}{1-\alpha }}{b}_{1}sin\left(\sqrt{\frac{1}{1-\alpha }}(x+A)\right),A\in \mathbb{R}.$

Solution 31.

If $\alpha <1$, $a_0=0$, $a_1=a_1$, $a_2=-\frac{2b_1}{\alpha -1}$, $b_0=0$, $b_1\ne 0$, $\lambda =\frac{a_1-\alpha a_1}{b_1}$, $\mu =-\frac{(\alpha -1)\left(-\alpha a_1^2+a_1^2+b_1^2\right)}{4b_1^2}$ and $\omega =0$, from which we obtain the solution

(46) $u_{31}\left(x\right)=\left(-(1-\alpha )a_{1}cos\left(\frac{1}{2}\sqrt{1-\alpha }(A+x)\right)+b_1\sqrt{1-\alpha }sin\left(\frac{1}{2}\sqrt{1-\alpha }(A+x)\right)\right){}^{\frac{2}{1-\alpha }},A\in \mathbb{R}.$

Solution 32.

If $\alpha <1$, $a_0\ne 0$, $a_1=0$, $a_2=-\frac{a_0}{\mu }$, $b_0=0$, $b_1=-\frac{a_0}{\mu }$, $\lambda =0$, $\mu =-\frac{1}{4(\alpha -1)}$ and $\omega =0$, from which we obtain the solution

(47) $u_{32}\left(x\right)=\frac{1}{2}sin\left(\sqrt{\frac{1}{1-\alpha }}(x+A)\right),A\in \mathbb{R}.$

Solution 33.

If $\alpha =\frac{1}{3}$, $a_0=0$, $a_1=a_1$, $a_2=3b_1$, $b_0=0$, $b_1\ne 0$, $\lambda =\frac{2a_1}{3b_1}$, $\mu =\frac{2a_1^2+3b_1^2}{18b_1^2}$ and $\omega =0$, from which we obtain the solution

(48) $u_{33}\left(x\right)={\left(2a_{1}cos\left(\frac{x+A}{\sqrt{6}}\right)-\sqrt{6}{b}_{1}sin\left(\frac{x+A}{\sqrt{6}}\right)\right)}^3,A\in \mathbb{R}.$

Solution 34.

If $\alpha =\frac{1}{2}$, $a_0\ne 0$, $a_1=0$, $a_2=-\frac{3a_0}{\mu }$, $b_0=0$, $b_1=-\frac{a_0}{\mu }$, $\lambda =0$, $\mu =\frac{1}{8}$ and $\omega =0$, from which we obtain the solution

(49) $u_{34}\left(x\right)={sin}^3\left(\frac{x+A}{2\sqrt{2}}\right)cos\left(\frac{x+A}{2\sqrt{2}}\right),A\in \mathbb{R}.$

Solution 35.

If $\alpha =\frac{1}{2}$, $a_0\ne 0$, $a_1=0$, $a_2=-\frac{a_0}{3\mu }$, $b_0=0$, $b_1=-\frac{a_0}{3\mu }$, $\lambda =0$, $\mu =\frac{1}{8}$ and $\omega =0$, from which we obtain the solution

(50) $u_{35}\left(x\right)=\frac{1}{8}\left(sin\left(\sqrt{2}(A+x)\right)+2sin\left(\frac{x+A}{\sqrt{2}}\right)\right),A\in \mathbb{R}.$

Three of the five solutions in the instability region are time-independent, and two are time-dependent. With these solutions, one of the difficulties is that they involve increasingly more parameters external to the model.

Solution 36.

If $\alpha >1$, $a_0=0$, $a_1=-\frac{b_1}{\sqrt{\alpha -1}}$, $a_2=-\frac{2b_1}{\alpha -1}$, $b_0=0$, $b_1\ne 0$, $\lambda =\sqrt{\alpha -1}$, $\mu =0$ and $\omega =0$, from which we obtain the solution

(51) $u_{36}\left(x\right)=e^{\frac{x}{\sqrt{\alpha -1}}}{\left(1-e^{\sqrt{\alpha -1}(A+x)}\right)}^{\frac{2}{1-\alpha }},A\in \mathbb{R}.$

Solution 37.

If $\alpha =\frac{a_1^2+b_1^2}{a_1^2}$, $a_0=0$, $a_1\ne 0$, $a_2=b_1$, $b_0=0$, $b_1\ne 0$, $\lambda =\frac{2a_1}{b_1}$, $\mu =0$ and $\omega =0$, from which we obtain the solution

(52) $u_{37}\left(x\right)=e^{-\frac{\lambda x}{2}}\left(1-e^{\lambda (x+A)}\right),A\in \mathbb{R}.$

Solution 38.

If $\alpha =\frac{a_1^2+b_1^2}{a_1^2}$, $a_0=0$, $a_1\ne 0$, $a_2=\frac{2b_1}{1-\alpha }$, $b_0=0$, $b_1\ne 0$, $\lambda =\frac{a_1-\alpha a_1}{b_1}$, $\mu =-\frac{(\alpha -1)\left(-\alpha a_1^2+a_1^2+b_1^2\right)}{4b_1^2}$ and $\omega =0$, from which we obtain the solution

(53) $u_{38}\left(x\right)=\frac{e^{\frac{a_1}{b_1}x}}{{\left(1-e^{\frac{a_1}{b_1}(x+A)}\right)}^{\frac{2a_1^2}{b_1^2}}},A\in \mathbb{R}.$

Solution 39.

If $\alpha =3$, we have $a_0=\frac{b_1}{6}$, $a_1=\mp \frac{b_1}{\sqrt{6}}$, $a_2=b_1$, $b_0=0$, $b_1\ne 0$, $\lambda =\mp \sqrt{\frac{2}{3}}$, $\mu =\frac{1}{6}$, and $\omega =\pm \frac{1}{3}\sqrt{\frac{2b_1}{3}}$, from which we obtain the solution

(54) $u_{39}(x,t)=\frac{3e^{\pm (\frac{A}{\sqrt{6}}-\frac{t}{9}+\frac{x}{\sqrt{6}})}\left(3\sqrt{6}A\mp 2t+3\sqrt{6}x\mp 18\right)}{9A\mp \sqrt{6}t+9x},A\in \mathbb{R}.$

Solution 40.

If $\alpha \ne 1$, $a_0=a_0$, $a_1\ne 0$, $b_0=0$, $b_1=0$, $a_2=b_2$, $\lambda =\frac{a_1}{\alpha -1}$, $\mu =\frac{(\alpha -1)}{2}$ and $\omega =a_0^3(\alpha -1)+a_0$, from which we obtain the solution

(55) $u_{40}(x,t)=A{e}^{\left(a_{0}x+a_0^{2}t(a_0^2(\alpha -1)-1)\right)},A\in \mathbb{R}.$

These solutions are combinations of exponentials; however, for several parameter values, no spike or wavelet patterns occur, and as a result, these combinations cannot be considered as a distribution.

4. Summary and Conclusions

In this work, we find several families of analytical solutions to the CWK equation for different values of the parameters, not only in the instability region. For this purpose, we use the generalized Kudryashov method, the $e^{-R\left(\xi \right)}$-expansion, and exponential function methods, which allows us to find exact solutions of nonlinear differential equations, including those with variable coefficients, non-integer powers, singular perturbation problems, and non-polynomial nonlinearities. These methods allow us to find analytical solutions to the CWK equation that have not been previously reported in the literature. Specifically, when $\alpha \ne 0$, the nonlinearity in the CWK equation cannot be analyzed with conventional methods. Our solutions are important because they provide insights into the behavior of the system and can be used in numerical simulations or experiments. Moreover, they can be used to develop new mathematical tools and techniques for analyzing and solving other nonlinear differential equations in fluid mechanics and related fields.

We found twenty-seven analytical solutions in the case $\alpha \ne 0$ by using the three methods. In the $\alpha <1$ case, the nonlinear term can be considered a perturbation of the linear behavior where inverse diffusion and fourth-order terms drive the dynamics. Although structure formation is not expected here, we find several wave-like stationary solutions. Thus, this region can be considered a pre-pattern formation region, and the solutions are budding patterns. The instability region $\alpha >1$ is where this model’s aggregate formation is expected. We first recover the stationary solution found in [14] and a pulse propagating without deformation by the generalized Kudryashov method. Some similar spike-like solutions were found by the $e^{-R}$-expansion method. The exponential function method did not yield any new physical solutions beyond those obtained by previous methods. This is not surprising, given that N.A. Kudryashov observed this in [27].

It is worth emphasizing that these methods allow us to obtain a wide range of behaviors, some previously obtained and that qualitatively resemble what was seen in the experiments. Furthermore, reverse diffusion and fourth-order terms still need to be explored since the traveling solutions were mainly found in this regime. The pattern formation and pulse motion mechanism could be better understood by studying each case in depth. On the other hand, the shortcoming of these methods is the large number of free parameters that appear in the solutions. One only way to fix them is to consider initial and boundary value problems, making the solutions more realistic and closer to the experimental phototaxis conditions, even in simple models like the one-dimensional CWK equation. We strongly believe the collection of physically meaningful solutions can guide the study of bacterial aggregate formation in phototaxis.

Author Contributions

Methodology, Software, Formal analysis, Writing—Original Draft, Visualization, A.L.-R.; Conceptualization, Methodology, Investigation, Supervision, O.G.-G.; Conceptualization, Methodology, Investigation, Supervision, G.C.-A. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Acknowledgments

ALR declares that he received support from the Mexican Council of Science and Technology (CONACyT), through scholarship No. 798462.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Solution 4 for [Forumla omitted. See PDF.]. Although no aggregate is produced, the wavefront propagates to the right with velocity [Forumla omitted. See PDF.].

Figure 2. Solution 8 reproduces the stationary finger-like distribution obtained in [14]. Values of [Forumla omitted. See PDF.] are presented. As [Forumla omitted. See PDF.] grows, the distribution becomes increasingly wider.

Figure 3. Solution 11 where a small pulse propagates to the right with velocity [Forumla omitted. See PDF.].

Figure 4. Solution 14 for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. Propagating wave front behavior is observed.

Figure 5. Solution 16 for fixed [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.], where increasing the value of [Forumla omitted. See PDF.] increases the amplitude of each curve.

Figure 6. Solution 18 for fixed [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. We show tree cases for [Forumla omitted. See PDF.], as [Forumla omitted. See PDF.] increases the amplitude of the spike decreases.

Appendix A. Stability Analysis

In [14], the authors derive Equation (1) from the continuum limit of the system of differential equations over a lattice formed by n bins, which is given by (A1)$$\left\{\begin{array}{c}\frac{d{R}_j}{dt}=a{R}_{j-1}-(a+c)R_j+c{U}_{j-1}{\eta }_{j-1}^{+}\hfill \\ \frac{d{L}_j}{dt}=a{L}_{j+1}-(a+c)L_j+c{U}_{j+1}{\eta }_{j+1}^{-};\hfill \end{array}\right.$$ where (A2)$${\eta }_j^{\pm }=\frac{{\sum }_{k=1}^d{U}_{j\pm k}}{{\sum }_{k=1}^d(U_{j+k}+U_{j-k})}and{U}_j=L_j+R_j.$$Here, $R_j\left(t\right)$ and $L_j\left(t\right)$ represent the density of right- and left-moving bacteria in the bin $j=1,...,n$ at time t; a represents the rate at which the bacterium moves one bin according to its orientation; and c represents the rate at which the bacterium moves after transitioning to a new orientation. The parameter d is the sensing radius of the bacterium.

Considering that $U_j=L_j+R_j$ and adding the two equalities of the system (A1), we obtain (A3)$$\frac{d{U}_j}{dt}=a(R_{j-1}+L_{j+1})-(a+c)U_j+c(U_{j-1}{\eta }_{j-1}^{+}+U_{j+1}{\eta }_{j+1}^{-}).$$A closed system can be constructed by defining (A4)$$V_j=R_{j+1}+L_{j-1}.$$Observe that immediately from system (A1) and (A2) we have (A5)$$\left\{\begin{array}{c}\frac{d{R}_{j+1}}{dt}=a{R}_j-(a+c)R_{j+1}+c{U}_j{\eta }_j^{+}\hfill \\ \frac{d{L}_{j-1}}{dt}=a{L}_j-(a+c)L_{j-1}+c{U}_j{\eta }_j^{-}.\hfill \end{array}\right.$$Now, adding the two equations in (A5), we obtain (A6)$$\frac{d{V}_j}{dt}=(a+c)(U_j-V_j).$$Now, rewriting (A7)$$R_{j-1}+L_{j+1}=U_{j-1}+U_{j+1}-V_j,$$ we conclude that the system (A1) is equivalent to (A8)$$\left\{\begin{array}{c}\frac{d{U}_j}{dt}=a(U_{j-1}+U_{j+1}-V_j)-(a+c)U_j+c(U_{j-1}{\eta }_{j-1}^{+}+U_{j+1}{\eta }_{j+1}^{-})\hfill \\ \frac{d{V}_j}{dt}=(a+c)(U_j-V_j).\hfill \end{array}\right.$$Model (A1) clearly allows for a homogeneous equilibrium $L_j=R_j=C$ for any constant C. Now, we will analyze the stability of this equilibrium. It is easier and more practical to conduct the analysis for the system (A8) whose steady state is given by $U_j=V_j=V$.

Consider the following perturbations:(A9)$$U_j=V+{\xi }_j\left(t\right);V_j=V+{\rho }_j\left(t\right),j=1,2,\dots ,n$$ where $|{\xi }_j|,|{\rho }_j|\ll 1$. We now obtain the linearized system from the system (A8):(A10)$$\begin{array}{c}\frac{d{\xi }_k}{dt}=(a+c/2)({\xi }_{k-1}+{\xi }_{j+1})-a{\rho }_j-(a+c){\xi }_j\hfill \\ \hfill +\frac{c}{2d}(2{\xi }_j+{\xi }_{j-1}+{\xi }_{j+1}-{\xi }_{j+d}-{\xi }_{j-d}-{\xi }_{j+d+1}-{\xi }_{j-d-1})\end{array}$$ (A11)$$\frac{d{\rho }_j}{dt}=(a+c)({\xi }_j-{\rho }_j).$$This $\left(2n\right)\times \left(2n\right)$ linear problem can be divided down into n subproblems of $2\times 2$. Make an ansatz (A12)$${\xi }_j=\xi e^{\lambda t}{e}^{\frac{2\pi mji}{n}};{\rho }_j=\rho e^{\lambda t}{e}^{\frac{2\pi mji}{n}},m=0,1,\dots ,n-1$$ to obtain (A13)$$\begin{array}{c}\lambda \xi =(2a+c)\xi cos\left(\theta \right)-a\rho -(a+c)\xi \hfill \\ \hfill +\frac{c}{2d}\xi \left(1+cos\left(\theta \right)-cos\left(d\theta \right)-cos\left(\right(d+1\left)\theta \right)\right)\end{array}$$ (A14)$$\lambda \rho =(a+c)(\xi -\rho )$$Here, $\theta =\frac{2\pi m}{n}$ with $m=0,1,\dots ,n-1$.

There are two eigenvalues for each possible value of m, for a total of $2n$ eigenvalues. The quadratic equation (A15)$${\lambda }^2-(g\left(\theta \right)-c)\lambda -(a+c)g\left(\theta \right)=0$$ gives the solution to the $2\times 2$ eigenvalue Problems (A13) and (A14), where function g is defined as (A16)$$g\left(\theta \right)=(2a+c)\left(cos\left(\theta \right)-1\right)+\frac{c}{2d}\left(1+cos\left(\theta \right)-cos\left(d\theta \right)-cos\left(\right(d+1\left)\theta \right)\right).$$Note that $g\left(\theta \right)-c\le 0$ for all $\theta$ so that a sufficient and necessary condition for stability is that $g\left(\theta \right)<0$ for all $\theta$.

Computations reveal that the instability occurs for the first time at $\theta =0$. Because $g\left(0\right)=g^{\prime }\left(0\right)=0$, the value of the threshold can be determined by setting $g^{″}\left(0\right)=dc-2a=0$. Consequently, we conclude that the critical value of the threshold is $c_0=\frac{2a}{d}$. The homogeneous steady state is therefore stable when $c<c_0$ and unstable when $c>c_0$, i.e., it is unstable if $c>\frac{2a}{d}$. The conclusion is obtained by spectral equivalence.

Appendix B. Separation of Variables Method for the Linear Case α = 0

Here, we discuss the case $\alpha =0$ occurring when, after changing orientation, the bacterium stops $c=0$, or when the rate of motion without changing orientation is very large $a\to \infty$, both for all finite d. In this case, Equation (1) reduces to the linear differential equation $u_t=-u_{xx}-u_{xxxx}$, i.e., only the reverse diffusion and long-range terms. The aggregate size is controlled by $\alpha$, so we cannot refer here to finger-like solutions. This equation can be solved by the well-known method of separation of variables where $u(x,t)=f\left(x\right)g\left(t\right)$ leads to the following:(A17)$$\frac{g^{\prime }}{g}=-\frac{(f^{″}+f^{″″})}{f}={\gamma }^2,$$ where the solutions can be directly obtained by considering the three possible cases for the separation parameter ${\gamma }^2$:

Case A1.

${\gamma }^2=0$ (A18) $$g^{\prime }=0and{f}^{″″}+f^{″}=0.$$ solving two linear differential equations gives (A19) $$g\left(t\right)=c_{1}andf\left(x\right)=k_1+k_{2}x+k_{3}cos\left(x\right)+k_{4}sin\left(x\right)$$ from which we obtain the family of solutions (A20) $$u_A(x,t)=C_1+C_{2}x+C_{3}cos\left(x\right)+C_{4}sin\left(x\right),C_1,C_2,C_3,C_4\in \mathbb{R}.$$

Case A2.

${\gamma }^2>0$ (A21) $$g^{\prime }-{\gamma }^{2}g=0and{f}^{″″}+f^{″}+{\gamma }^{2}f=0.$$ solving two linear differential equations gives (A22) $$g\left(t\right)=c_1{e}^{{\gamma }^{2}t}andf\left(x\right)=k_1{e}^{\left(\frac{-1-\beta }{2}\right)x}+k_2{e}^{\left(\frac{-1+\beta }{2}\right)x}+k_{3}x{e}^{\left(\frac{-1-\beta }{2}\right)x}+k_{4}x{e}^{\left(\frac{-1+\beta }{2}\right)x}$$ from which, considering $\beta =\sqrt{1-4{\gamma }^2}$, we obtain the family of solutions (A23) $$u_B(x,t)=e^{{\gamma }^{2}t}\left(C_1{e}^{\left(\frac{-1-\beta }{2}\right)x}+C_2{e}^{\left(\frac{-1+\beta }{2}\right)x}+C_{3}x{e}^{\left(\frac{-1-\beta }{2}\right)x}+C_{4}x{e}^{\left(\frac{-1+\beta }{2}\right)x}\right),C_1,C_2,C_3,C_4\in \mathbb{R}.$$

Case A2.

${\gamma }^2<0$ (A24) $$g^{\prime }+{\gamma }^{2}g=0and{f}^{″″}+f^{″}-{\gamma }^{2}f=0.$$ solving two linear differential equations gives (A25) $$g\left(t\right)=c_1{e}^{-{\gamma }^{2}t}andf\left(x\right)=k_1{e}^{\left(\frac{-1-\delta }{2}\right)x}+k_2{e}^{\left(\frac{-1+\delta }{2}\right)x}+k_{3}x{e}^{\left(\frac{-1-\delta }{2}\right)x}+k_{4}x{e}^{\left(\frac{-1+\delta }{2}\right)x}$$ from which, considering $\delta =\sqrt{1+4{\gamma }^2}$, we obtain the family of solutions (A26) $$u_C(x,t)=e^{-{\gamma }^{2}t}\left(C_1{e}^{\left(\frac{-1-\delta }{2}\right)x}+C_2{e}^{\left(\frac{-1+\delta }{2}\right)x}+C_{3}x{e}^{\left(\frac{-1-\delta }{2}\right)x}+C_{4}x{e}^{\left(\frac{-1+\delta }{2}\right)x}\right),C_1,C_2,C_3,C_4\in \mathbb{R}.$$

All five constants in each case can be fixed through the corresponding initial and boundary conditions. We generally observe that, according to the sign of ${\gamma }^2$, we can have oscillatory solutions or increasing and decreasing exponential solutions. Consider also that the principle of superposition of solutions is valid for the present linear case. We have one family of standing wave-like solutions and two families of traveling wavefront-like solutions, similar to those obtained with the previous methods. Moreover, these solutions coincide with those obtained with the proposed methods. We show some examples of this below.

Example A1.

If we consider the family of solutions $u_1$: $$u_1(x,t)=cosh(2t-x)-sinh(2t-x)=e^{-(2t-x)}+1,$$ this is derived from the present method using $u_C$ and considering $C_1=0$, $C_2=1$, $\delta -1=2$, ${\gamma }^2=4$, $C_3=C_4=0$.

Example A2.

If we consider the family of solutions $u_3$: $$u_3(x,t)=cosh(20t+2x)-sinh(20t+2x)=e^{-(2x+20t)}-1,$$ this is derived from the present method using $u_C$ and considering $C_1=0$, $C_2=1$, $\delta -1=-4$, ${\gamma }^2=20$, $C_3=C_4=0$.

Example A3.

If we consider the family of solutions $u_{13}$: $$u_{13}\left(x\right)=e^{-\lambda \left(({\lambda }^3+\lambda )t+x\right)}-e^{\lambda A},$$ this is derived from the present method using $u_C$ and considering $C_1=1$, $C_2=0$, $\delta +1=2\lambda$, ${\gamma }^2={\lambda }^4+{\lambda }^2$, $C_3=C_4=0$; moreover, $e^{\lambda A}$ is a constant.

Example A4.

If we consider the family of solutions $u_{28}$: $$u_{28}(x,t)=sinh\left(-\lambda A+{\lambda }^2({\lambda }^2+1)t-\lambda x\right)-cosh\left(-\lambda A+{\lambda }^2({\lambda }^2+1)t-\lambda x\right)+1,$$ this is derived from the present method using $u_C$ and considering $C_1=0$, $C_2=-e^{\lambda A}$, $\delta -1=2\lambda$, ${\gamma }^2={\lambda }^2({\lambda }^2+1)$, $C_3=C_4=0$.

The other families for $\alpha =0$ are obtained similarly and consider the superposition principle. In the present case of $\alpha =0$, aggregate formation is not expected when the nonlinear term does not appear in the CWK equation. However, among the solutions found are time-propagating wavefront solutions and some standing wave solutions that could be interpreted as finger-shaped distributions. Such solutions are interesting but cannot be considered a final steady state; they must start with that form.

Unfortunately, for $\alpha \ne 0$, Equation (1) is nonlinear. Note that when the derivative is expanded, the nonlinearity becomes more involved:$$u_t=-u_{xx}-u_{xxxx}+\alpha \left[\frac{u{u}_x{u}_{xxx}+u{u}_{xx}^2-u_x^2{u}_{xx}}{u^2}\right].$$

This equation cannot always be solved by the method of separation of variables, nor can all the families of solutions found in this work be obtained.

Appendix C. Algebraic System for Kudryashov Method

The Kudryashov method algorithm requires solving the following system of equations for the unknowns $\{a_0,a_1,a_2,b_0,b_1,\omega \}$. $$\begin{array}{c}\hfill Q^0:\alpha a_0^4+a_0{b}_0^3\omega -a_0^2{b}_0^2-a_0^4=0,\end{array}$$ $$\begin{array}{cc}\hfill Q^1:& 4\alpha a_1{a}_0^3+4\alpha a_0^3{b}_1-4\alpha a_1{a}_0^2{b}_0-\alpha a_0^2{b}_0{b}_1+\alpha a_1{a}_0{b}_0^2+3a_0{b}_0^2{b}_1\omega +a_1{b}_0^3\omega -6a_0^3{b}_1\hfill \\ \hfill & +6a_1{a}_0^2{b}_0+2a_0^2{b}_0{b}_1-6a_1{a}_0{b}_0^2-2a_0{b}_0^2{b}_1+2a_1{b}_0^3-4a_1{a}_0^3=0,\hfill \end{array}$$ (A27)$$\begin{array}{cc}\hfill Q^2:& 4\alpha a_2{a}_0^3+6\alpha a_1^2{a}_0^2-4\alpha a_0^3{b}_1+2\alpha a_0^2{b}_1^2+4\alpha a_1{a}_0^2{b}_0-8\alpha a_2{a}_0^2{b}_0+8\alpha a_1{a}_0^2{b}_1+3\alpha a_0^2{b}_0{b}_1\hfill \\ \hfill & -3\alpha a_1{a}_0{b}_0^2+4\alpha a_2{a}_0{b}_0^2-8\alpha a_1^2{a}_0{b}_0-4\alpha a_1{a}_0{b}_0{b}_1+2\alpha a_1^2{b}_0^2+3a_0{b}_0{b}_1^2\omega \hfill \\ \hfill & +a_2{b}_0^3\omega +3a_1{b}_0^2{b}_1\omega +6a_0^3{b}_1-8a_0^2{b}_1^2-6a_1{a}_0^2{b}_0+12a_2{a}_0^2{b}_0-12a_1{a}_0^2{b}_1\hfill \\ \hfill & -12a_0^2{b}_0{b}_1+12a_1{a}_0{b}_0^2-18a_2{a}_0{b}_0^2+2a_0{b}_0{b}_1^2+12a_1^2{a}_0{b}_0+8a_0{b}_0^2{b}_1\hfill \\ \hfill & +10a_1{a}_0{b}_0{b}_1-8a_1{b}_0^3+10a_2{b}_0^3-8a_1^2{b}_0^2-2a_1{b}_0^2{b}_1-4a_2{a}_0^3-6a_1^2{a}_0^2=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill Q^3:& 4\alpha a_0{a}_1^3+12\alpha a_0^2{a}_2{a}_1-4\alpha a_1^3{b}_0-5\alpha a_1^2{b}_0^2+8\alpha a_0{a}_1^2{b}_0+4\alpha a_0{a}_1^2{b}_1-\alpha a_1^2{b}_0{b}_1\hfill \\ \hfill & +2\alpha a_0{a}_1{b}_0^2+9\alpha a_2{a}_1{b}_0^2+\alpha a_0{a}_1{b}_1^2-24\alpha a_0{a}_2{a}_1{b}_0-8\alpha a_0^2{a}_1{b}_1+8\alpha a_0{a}_1{b}_0{b}_1\hfill \\ \hfill & -10\alpha a_0{a}_2{b}_0^2-3\alpha a_0^2{b}_1^2+8\alpha a_0^2{a}_2{b}_0+4\alpha a_0^2{a}_2{b}_1-2\alpha a_0^2{b}_0{b}_1-2\alpha a_0{a}_2{b}_0{b}_1\hfill \\ \hfill & +3a_1{b}_0{b}_1^2\omega +a_0{b}_1^3\omega +3a_2{b}_0^2{b}_1\omega +6a_1^3{b}_0+18a_1^2{b}_0^2-12a_0{a}_1^2{b}_0-6a_0{a}_1^2{b}_1\hfill \\ \hfill & +2a_1^2{b}_0{b}_1+12a_1{b}_0^3-8a_0{a}_1{b}_0^2-34a_2{a}_1{b}_0^2-6a_0{a}_1{b}_1^2+2a_1{b}_0{b}_1^2+36a_0{a}_2{a}_1{b}_0\hfill \\ \hfill & +12a_0^2{a}_1{b}_1+8a_1{b}_0^2{b}_1-28a_0{a}_1{b}_0{b}_1-40a_2{b}_0^3-2a_0{b}_1^3+40a_0{a}_2{b}_0^2+10a_0^2{b}_1^2-8a_0{b}_0{b}_1^2\hfill \\ \hfill & -12a_0^2{a}_2{b}_0-12a_0{b}_0^2{b}_1+10a_2{b}_0^2{b}_1-6a_0^2{a}_2{b}_1+8a_0^2{b}_0{b}_1-4a_0{a}_1^3-12a_0^2{a}_2{a}_1=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill Q^4:& \alpha a_1^4+12\alpha a_0{a}_2{a}_1^2+6\alpha a_0^2{a}_2^2+4\alpha a_1^3{b}_0+3\alpha a_1^2{b}_0^2-16\alpha a_2{a}_1^2{b}_0-4\alpha a_0{a}_1^2{b}_1\hfill \\ \hfill & +\alpha a_1^2{b}_0{b}_1-21\alpha a_2{a}_1{b}_0^2-\alpha a_0{a}_1{b}_1^2+24\alpha a_0{a}_2{a}_1{b}_0-4\alpha a_0{a}_1{b}_0{b}_1-18a_2{a}_1{b}_0{b}_1\hfill \\ \hfill & -4a_1^2{b}_0{b}_1+4\alpha a_2{a}_1{b}_0{b}_1+8\alpha a_2^2{b}_0^2+6\alpha a_0{a}_2{b}_0^2+\alpha a_0^2{b}_1^2-16\alpha a_0{a}_2^2{b}_0-4\alpha a_0^2{a}_2{b}_1\hfill \\ \hfill & +2\alpha a_0{a}_2{b}_0{b}_1+a_1{b}_1^3\omega +3a_2{b}_0{b}_1^2\omega -6a_1^3{b}_0-11a_1^2{b}_0^2-a_1^2{b}_1^2+24a_2{a}_1^2{b}_0+6a_0{a}_1^2{b}_1\hfill \\ \hfill & -6a_1{b}_0^3+76a_2{a}_1{b}_0^2+4a_0{a}_1{b}_1^2-2a_1{b}_0{b}_1^2-36a_0{a}_2{a}_1{b}_0-6a_1{b}_0^2{b}_1+14a_0{a}_1{b}_0{b}_1\hfill \\ \hfill & +54a_2{b}_0^3+2a_0{b}_1^3-29a_2^2{b}_0^2-24a_0{a}_2{b}_0^2-3a_0^2{b}_1^2-4a_0{a}_2{b}_1^2+6a_0{b}_0{b}_1^2+8a_2{b}_0{b}_1^2\hfill \\ \hfill & +24a_0{a}_2^2{b}_0+6a_0{b}_0^2{b}_1-46a_2{b}_0^2{b}_1+6a_0^2{a}_2{b}_1-a_1^4-12a_0{a}_2{a}_1^2-6a_0^2{a}_2^2=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill Q^5:& 4\alpha a_2{a}_1^3+12\alpha a_0{a}_2^2{a}_1+16\alpha a_2{a}_1^2{b}_0-4\alpha a_2{a}_1^2{b}_1+12\alpha a_2{a}_1{b}_0^2+\alpha a_2{a}_1{b}_1^2\hfill \\ \hfill & +44a_2{a}_1{b}_0{b}_1-20\alpha a_2^2{a}_1{b}_0-12\alpha a_2{a}_1{b}_0{b}_1-18\alpha a_2^2{b}_0^2+16\alpha a_0{a}_2^2{b}_0-4\alpha a_0{a}_2^2{b}_1\hfill \\ \hfill & +7\alpha a_2^2{b}_0{b}_1+a_2{b}_1^3\omega -24a_2{a}_1^2{b}_0+6a_2{a}_1^2{b}_1-44a_2{a}_1{b}_0^2-6a_2{a}_1{b}_1^2+30a_2^2{a}_1{b}_0\hfill \\ \hfill & -24a_2{b}_0^3+2a_2{b}_1^3+64a_2^2{b}_0^2+4a_0{a}_2{b}_1^2-32a_2{b}_0{b}_1^2-24a_0{a}_2^2{b}_0+6a_0{a}_2^2{b}_1\hfill \\ \hfill & +72a_2{b}_0^2{b}_1-26a_2^2{b}_0{b}_1-4a_0{a}_2{b}_0{b}_1-4a_2{a}_1^3-12a_0{a}_2^2{a}_1=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill Q^6:& 4\alpha a_0{a}_2^3+6\alpha a_1^2{a}_2^2-8\alpha a_2^3{b}_0+10\alpha a_2^2{b}_0^2+2\alpha a_2^2{b}_1^2+20\alpha a_1{a}_2^2{b}_0+12a_2^3{b}_0\hfill \\ \hfill & +4\alpha a_0{a}_2^2{b}_1-8\alpha a_1{a}_2^2{b}_1-17\alpha a_2^2{b}_0{b}_1-3\alpha a_1{a}_2{b}_1^2+4\alpha a_1^2{a}_2{b}_1+8\alpha a_1{a}_2{b}_0{b}_1\hfill \\ \hfill & -36a_2^2{b}_0^2-8a_2^2{b}_1^2-30a_1{a}_2^2{b}_0-6a_0{a}_2^2{b}_1+12a_1{a}_2^2{b}_1+60a_2^2{b}_0{b}_1-8a_2{b}_1^3-2a_0{a}_2{b}_1^2\hfill \\ \hfill & +12a_1{a}_2{b}_1^2+48a_2{b}_0{b}_1^2-6a_1^2{a}_2{b}_1-36a_2{b}_0^2{b}_1-30a_1{a}_2{b}_0{b}_1-4a_0{a}_2^3-6a_1^2{a}_2^2=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill Q^7:& 4\alpha a_1{a}_2^3+8\alpha a_2^3{b}_0-4\alpha a_2^3{b}_1-5\alpha a_2^2{b}_1^2+8\alpha a_1{a}_2^2{b}_1+10\alpha a_2^2{b}_0{b}_1+2\alpha a_1{a}_2{b}_1^2-12a_2^3{b}_0\hfill \\ \hfill & +6a_2^3{b}_1+18a_2^2{b}_1^2-12a_1{a}_2^2{b}_1-36a_2^2{b}_0{b}_1+12a_2{b}_1^3-8a_1{a}_2{b}_1^2-24a_2{b}_0{b}_1^2-4a_1{a}_2^3=0,\hfill \end{array}$$ $$\begin{array}{c}\hfill Q^8:\alpha a_2^4+4\alpha a_2^3{b}_1+3\alpha a_2^2{b}_1^2-6a_2^3{b}_1-11a_2^2{b}_1^2-6a_2{b}_1^3-a_2^4=0.\end{array}$$

Appendix D. Algebraic System for the e^−R(ξ) -Expansion Method

In the algorithm of the $e^{-R\left(\xi \right)}$-expansion method, the algorithm needs to solve the following system of equations for the $a_0,a_1$, and $\omega$; here we present the system for each power of $e^{-R\left(\xi \right)}$ from 0 to 4. $$\begin{array}{cc}\hfill e^0:& a_1\mu \left((\alpha -3)a_1\mu +{\lambda }^2+2\mu +1\right)+a_0\left((\alpha -4)a_1\lambda \mu +\omega \right)\hfill \\ \hfill & +a_0^2\left((6-4\alpha )a_1\mu -1\right)+(\alpha -1)a_0^4=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill e^{-R\left(z\right)}:& \alpha a_0{a}_1{\lambda }^2+3\alpha a_1^2\lambda \mu -4\alpha a_0^2{a}_1\lambda -8\alpha a_0{a}_1^2\mu +2\alpha a_0{a}_1\mu \hfill \\ \hfill & +4\alpha a_0^3{a}_1+a_1{\lambda }^3-4a_0{a}_1{\lambda }^2-10a_1^2\lambda \mu +8a_1\lambda \mu +6a_0^2{a}_1\lambda \hfill \\ \hfill & +a_1\lambda +12a_0{a}_1^2\mu -8a_0{a}_1\mu +a_1\omega -4a_0^3{a}_1-2a_0{a}_1=0,\hfill \end{array}$$ (A28)$$\begin{array}{cc}\hfill e^{-2R\left(z\right)}:& 2\alpha a_1^2{\lambda }^2-8\alpha a_0{a}_1^2\lambda +3\alpha a_0{a}_1\lambda -4\alpha a_1^3\mu +4\alpha a_1^2\mu \hfill \\ \hfill & +6\alpha a_0^2{a}_1^2-4\alpha a_0^2{a}_1-7a_1^2{\lambda }^2+7a_1{\lambda }^2+12a_0{a}_1^2\lambda -12a_0{a}_1\lambda \hfill \\ \hfill & +6a_1^3\mu -14a_1^2\mu +8a_1\mu -6a_0^2{a}_1^2-a_1^2+6a_0^2{a}_1+a_1=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill e^{-3R\left(z\right)}:& -4\alpha a_1^3\lambda +5\alpha a_1^2\lambda +4\alpha a_0{a}_1^3-8\alpha a_0{a}_1^2+2\alpha a_0{a}_1\hfill \\ \hfill & +6a_1^3\lambda -18a_1^2\lambda +12a_1\lambda -4a_0{a}_1^3+12a_0{a}_1^2-8a_0{a}_1=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill e^{-4R\left(z\right)}:& \alpha a_1^4-4\alpha a_1^3+3\alpha a_1^2-a_1^4+6a_1^3-11a_1^2+6a_1=0.\hfill \end{array}$$

Appendix E. System of Equations for Exponential Function Method

To obtain solution through the exponential function method, the following nonlinear algebraic system needs to be solved:$$\begin{array}{c}\hfill e^0:a_2\left(-4a_2{b}_1+a_2^2+3b_1^2\right)\left((\alpha -1)a_2+2b_1\right)=0,\end{array}$$ $$\begin{array}{cc}\hfill e^{Q\left(z\right)}:& 4\alpha a_1{a}_2^3-4\alpha a_2^3{b}_1\lambda +5\alpha a_2^2{b}_1^2\lambda -8\alpha a_2^3{b}_0-8\alpha a_1{a}_2^2{b}_1+10\alpha a_2^2{b}_0{b}_1+2\alpha a_1{a}_2{b}_1^2+6a_2^3{b}_1\lambda \hfill \\ \hfill & -18a_2^2{b}_1^2\lambda +12a_2{b}_1^3\lambda +12a_2^3{b}_0+12a_1{a}_2^2{b}_1-36a_2^2{b}_0{b}_1-8a_1{a}_2{b}_1^2+24a_2{b}_0{b}_1^2-4a_1{a}_2^3=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill e^{2Q\left(z\right)}:& 4\alpha a_0{a}_2^3+6\alpha a_1^2{a}_2^2+2\alpha a_2^2{b}_1^2{\lambda }^2-8\alpha a_2^3{b}_0\lambda -8\alpha a_1{a}_2^2{b}_1\lambda +17\alpha a_2^2{b}_0{b}_1\lambda +3\alpha a_1{a}_2{b}_1^2\lambda \hfill \\ \hfill & -4\alpha a_2^3{b}_1\mu +4\alpha a_2^2{b}_1^2\mu +10\alpha a_2^2{b}_0^2-20\alpha a_1{a}_2^2{b}_0-4\alpha a_0{a}_2^2{b}_1-4\alpha a_1^2{a}_2{b}_1+8\alpha a_1{a}_2{b}_0{b}_1\hfill \\ \hfill & -7a_2^2{b}_1^2{\lambda }^2+7a_2{b}_1^3{\lambda }^2+12a_2^3{b}_0\lambda +12a_1{a}_2^2{b}_1\lambda -60a_2^2{b}_0{b}_1\lambda -12a_1{a}_2{b}_1^2\lambda +48a_2{b}_0{b}_1^2\lambda \hfill \\ \hfill & +6a_2^3{b}_1\mu -14a_2^2{b}_1^2\mu +8a_2{b}_1^3\mu -36a_2^2{b}_0^2-a_2^2{b}_1^2+30a_1{a}_2^2{b}_0+6a_0{a}_2^2{b}_1+a_2{b}_1^3-2a_0{a}_2{b}_1^2\hfill \\ \hfill & +6a_1^2{a}_2{b}_1+36a_2{b}_0^2{b}_1-30a_1{a}_2{b}_0{b}_1-4a_0{a}_2^3-6a_1^2{a}_2^2=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill e^{3Q\left(z\right)}:& 4\alpha a_2{a}_1^3+12\alpha a_0{a}_2^2{a}_1+\alpha a_2{a}_1{b}_1^2{\lambda }^2+7\alpha a_2^2{b}_0{b}_1{\lambda }^2+3\alpha a_2^2{b}_1^2\lambda \mu -4\alpha a_2{a}_1^2{b}_1\lambda \hfill \\ \hfill & -20\alpha a_2^2{a}_1{b}_0\lambda +12\alpha a_2{a}_1{b}_0{b}_1\lambda +18\alpha a_2^2{b}_0^2\lambda -4\alpha a_0{a}_2^2{b}_1\lambda +2\alpha a_2{a}_1{b}_1^2\mu \hfill \\ \hfill & -8\alpha a_2^2{a}_1{b}_1\mu -8\alpha a_2^3{b}_0\mu +14\alpha a_2^2{b}_0{b}_1\mu -16\alpha a_2{a}_1^2{b}_0+12\alpha a_2{a}_1{b}_0^2-16\alpha a_0{a}_2^2{b}_0\hfill \\ \hfill & +a_2{b}_1^3{\lambda }^3-4a_2{a}_1{b}_1^2{\lambda }^2+28a_2{b}_0{b}_1^2{\lambda }^2-24a_2^2{b}_0{b}_1{\lambda }^2+8a_2{b}_1^3\lambda \mu -10a_2^2{b}_1^2\lambda \mu \hfill \\ \hfill & +6a_2{a}_1^2{b}_1\lambda +30a_2^2{a}_1{b}_0\lambda -44a_2{a}_1{b}_0{b}_1\lambda +a_2{b}_1^3\lambda -64a_2^2{b}_0^2\lambda -4a_0{a}_2{b}_1^2\lambda +6a_0{a}_2^2{b}_1\lambda \hfill \\ \hfill & +72a_2{b}_0^2{b}_1\lambda -8a_2{a}_1{b}_1^2\mu +12a_2^2{a}_1{b}_1\mu +32a_2{b}_0{b}_1^2\mu +12a_2^3{b}_0\mu -48a_2^2{b}_0{b}_1\mu \hfill \\ \hfill & +a_2{b}_1^3\omega +24a_2{a}_1^2{b}_0-44a_2{a}_1{b}_0^2-2a_2{a}_1{b}_1^2+24a_2{b}_0^3+4a_2{b}_0{b}_1^2\hfill \\ \hfill & +24a_0{a}_2^2{b}_0-2a_2^2{b}_0{b}_1-4a_0{a}_2{b}_0{b}_1-4a_2{a}_1^3-12a_0{a}_2^2{a}_1=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill e^{4Q\left(z\right)}:& \alpha a_1^4-a_1^4-4\alpha b_0{a}_1^3+6b_0{a}_1^3+3\alpha b_0^2{a}_1^2-11b_0^2{a}_1^2-b_1^2{a}_1^2+12\alpha a_0{a}_2{a}_1^2-12a_0{a}_2{a}_1^2\hfill \\ \hfill & -16\alpha \lambda a_2{b}_0{a}_1^2+24\lambda a_2{b}_0{a}_1^2+4\alpha a_0{b}_1{a}_1^2-6a_0{b}_1{a}_1^2-4\alpha \mu a_2{b}_1{a}_1^2+6\mu a_2{b}_1{a}_1^2\hfill \\ \hfill & -\alpha \lambda b_0{b}_1{a}_1^2+4\lambda b_0{b}_1{a}_1^2+6b_0^3{a}_1+\omega b_1^3{a}_1+21\alpha \lambda a_2{b}_0^2{a}_1-76\lambda a_2{b}_0^2{a}_1\hfill \\ \hfill & +\alpha \lambda a_0{b}_1^2{a}_1-4\lambda a_0{b}_1^2{a}_1+\alpha \lambda \mu a_2{b}_1^2{a}_1-4\lambda \mu a_2{b}_1^2{a}_1+{\lambda }^2{b}_0{b}_1^2{a}_1+2\mu b_0{b}_1^2{a}_1\hfill \\ \hfill & +b_0{b}_1^2{a}_1-20\alpha \mu a_2^2{b}_0{a}_1+30\mu a_2^2{b}_0{a}_1-24\alpha a_0{a}_2{b}_0{a}_1+36a_0{a}_2{b}_0{a}_1-6\lambda b_0^2{b}_1{a}_1\hfill \\ \hfill & -4\alpha a_0{b}_0{b}_1{a}_1+14a_0{b}_0{b}_1{a}_1+4\alpha {\lambda }^2{a}_2{b}_0{b}_1{a}_1-14{\lambda }^2{a}_2{b}_0{b}_1{a}_1+8\alpha \mu a_2{b}_0{b}_1{a}_1\hfill \\ \hfill & -28\mu a_2{b}_0{b}_1{a}_1-4a_2{b}_0{b}_1{a}_1+54\lambda a_2{b}_0^3-{\lambda }^2{a}_0{b}_1^3-2\mu a_0{b}_1^3-a_0{b}_1^3+2{\mu }^2{a}_2{b}_1^3\hfill \\ \hfill & +{\lambda }^2\mu a_2{b}_1^3+\mu a_2{b}_1^3+6\alpha a_0^2{a}_2^2-6a_0^2{a}_2^2+8\alpha {\lambda }^2{a}_2^2{b}_0^2-28{\lambda }^2{a}_2^2{b}_0^2+16\alpha \mu a_2^2{b}_0^2\hfill \\ \hfill & -56\mu a_2^2{b}_0^2-a_2^2{b}_0^2+6\alpha a_0{a}_2{b}_0^2-24a_0{a}_2{b}_0^2+\alpha a_0^2{b}_1^2-3a_0^2{b}_1^2+\alpha {\mu }^2{a}_2^2{b}_1^2\hfill \\ \hfill & -3{\mu }^2{a}_2^2{b}_1^2-2{\lambda }^2{a}_0{a}_2{b}_1^2-4\mu a_0{a}_2{b}_1^2-2a_0{a}_2{b}_1^2+6\lambda a_0{b}_0{b}_1^2+4{\lambda }^3{a}_2{b}_0{b}_1^2\hfill \\ \hfill & +4\lambda a_2{b}_0{b}_1^2+32\lambda \mu a_2{b}_0{b}_1^2+3\omega a_2{b}_0{b}_1^2-16\alpha \lambda a_0{a}_2^2{b}_0+24\lambda a_0{a}_2^2{b}_0\hfill \\ \hfill & -4\alpha \mu a_0{a}_2^2{b}_1+6\mu a_0{a}_2^2{b}_1-6a_0{b}_0^2{b}_1+41{\lambda }^2{a}_2{b}_0^2{b}_1+46\mu a_2{b}_0^2{b}_1+5a_2{b}_0^2{b}_1\hfill \\ \hfill & +4\alpha a_0^2{a}_2{b}_1-6a_0^2{a}_2{b}_1+11\alpha \lambda \mu a_2^2{b}_0{b}_1-36\lambda \mu a_2^2{b}_0{b}_1-2\alpha \lambda a_0{a}_2{b}_0{b}_1=0,\hfill \end{array}$$ (A29)$$\begin{array}{cc}\hfill e^{5Q\left(z\right)}:& -a_0{b}_1^3{\lambda }^3+a_1{b}_0{b}_1^2{\lambda }^3+5a_2{b}_0^2{b}_1{\lambda }^3+38a_2{b}_0^3{\lambda }^2+9\alpha a_1{a}_2{b}_0^2{\lambda }^2-32a_1{a}_2{b}_0^2{\lambda }^2\hfill \\ \hfill & +\alpha a_0{a}_1{b}_1^2{\lambda }^2-4a_0{a}_1{b}_1^2{\lambda }^2+10a_0{b}_0{b}_1^2{\lambda }^2+4\mu a_2{b}_0{b}_1^2{\lambda }^2-10a_1{b}_0^2{b}_1{\lambda }^2\hfill \\ \hfill & -\alpha a_1^2{b}_0{b}_1{\lambda }^2+4a_1^2{b}_0{b}_1{\lambda }^2-2\alpha a_0{a}_2{b}_0{b}_1{\lambda }^2+4a_0{a}_2{b}_0{b}_1{\lambda }^2+12a_1{b}_0^3\lambda \hfill \\ \hfill & -8\mu a_0{b}_1^3\lambda -a_0{b}_1^3\lambda +5\alpha a_1^2{b}_0^2\lambda -18a_1^2{b}_0^2\lambda +14\alpha \mu a_2^2{b}_0^2\lambda -48\mu a_2^2{b}_0^2\lambda \hfill \\ \hfill & +10\alpha a_0{a}_2{b}_0^2\lambda -40a_0{a}_2{b}_0^2\lambda +3\alpha a_0^2{b}_1^2\lambda -10a_0^2{b}_1^2\lambda -4\mu a_0{a}_2{b}_1^2\lambda +8\mu a_1{b}_0{b}_1^2\lambda \hfill \\ \hfill & +a_1{b}_0{b}_1^2\lambda -4\alpha a_1^3{b}_0\lambda +6a_1^3{b}_0\lambda -24\alpha a_0{a}_1{a}_2{b}_0\lambda +36a_0{a}_1{a}_2{b}_0\lambda +4\alpha a_0{a}_1^2{b}_1\lambda \hfill \\ \hfill & -6a_0{a}_1^2{b}_1\lambda -12a_0{b}_0^2{b}_1\lambda +40\mu a_2{b}_0^2{b}_1\lambda +5a_2{b}_0^2{b}_1\lambda +4\alpha a_0^2{a}_2{b}_1\lambda -6a_0^2{a}_2{b}_1\lambda \hfill \\ \hfill & -8\alpha a_0{a}_1{b}_0{b}_1\lambda +28a_0{a}_1{b}_0{b}_1\lambda +4\alpha \mu a_1{a}_2{b}_0{b}_1\lambda -12\mu a_1{a}_2{b}_0{b}_1\lambda +4\alpha a_0{a}_1^3\hfill \\ \hfill & -4a_0{a}_1^3+40\mu a_2{b}_0^3+2a_2{b}_0^3+\omega a_0{b}_1^3+2\alpha a_0{a}_1{b}_0^2-8a_0{a}_1{b}_0^2+18\alpha \mu a_1{a}_2{b}_0^2\hfill \\ \hfill & -64\mu a_1{a}_2{b}_0^2-2a_1{a}_2{b}_0^2+2\alpha \mu a_0{a}_1{b}_1^2-8\mu a_0{a}_1{b}_1^2-2a_0{a}_1{b}_1^2+8\mu a_0{b}_0{b}_1^2-2a_0{b}_0{b}_1^2\hfill \\ \hfill & +3\omega a_1{b}_0{b}_1^2+8{\mu }^2{a}_2{b}_0{b}_1^2+4\mu a_2{b}_0{b}_1^2+12\alpha a_0^2{a}_1{a}_2-12a_0^2{a}_1{a}_2-8\alpha a_0{a}_1^2{b}_0\hfill \\ \hfill & +12a_0{a}_1^2{b}_0-16\alpha \mu a_0{a}_2^2{b}_0+24\mu a_0{a}_2^2{b}_0-8\alpha a_0^2{a}_2{b}_0+12a_0^2{a}_2{b}_0-16\alpha \mu a_1^2{a}_2{b}_0\hfill \\ \hfill & +24\mu a_1^2{a}_2{b}_0-8\mu a_1{b}_0^2{b}_1+2a_1{b}_0^2{b}_1+3\omega a_2{b}_0^2{b}_1+8\alpha a_0^2{a}_1{b}_1-12a_0^2{a}_1{b}_1\hfill \\ \hfill & -2\alpha a_0^2{b}_0{b}_1+8a_0^2{b}_0{b}_1-2\alpha \mu a_1^2{b}_0{b}_1+8\mu a_1^2{b}_0{b}_1-2a_1^2{b}_0{b}_1+4\alpha {\mu }^2{a}_2^2{b}_0{b}_1\hfill \\ \hfill & -12{\mu }^2{a}_2^2{b}_0{b}_1-4\alpha \mu a_0{a}_2{b}_0{b}_1+8\mu a_0{a}_2{b}_0{b}_1-4a_0{a}_2{b}_0{b}_1=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill e^{6Q\left(z\right)}:& 8a_2{b}_0^3{\lambda }^3+4a_0{b}_0{b}_1^2{\lambda }^3-4a_1{b}_0^2{b}_1{\lambda }^3+7a_1{b}_0^3{\lambda }^2-7\mu a_0{b}_1^3{\lambda }^2+2\alpha a_1^2{b}_0^2{\lambda }^2-7a_1^2{b}_0^2{\lambda }^2\hfill \\ \hfill & +4\alpha a_0{a}_2{b}_0^2{\lambda }^2-16a_0{a}_2{b}_0^2{\lambda }^2+2\alpha a_0^2{b}_1^2{\lambda }^2-7a_0^2{b}_1^2{\lambda }^2+7\mu a_1{b}_0{b}_1^2{\lambda }^2-7a_0{b}_0^2{b}_1{\lambda }^2\hfill \\ \hfill & -\mu a_2{b}_0^2{b}_1{\lambda }^2-4\alpha a_0{a}_1{b}_0{b}_1{\lambda }^2+14a_0{a}_1{b}_0{b}_1{\lambda }^2+52\mu a_2{b}_0^3\lambda +2a_2{b}_0^3\lambda +3\alpha a_0{a}_1{b}_0^2\lambda \hfill \\ \hfill & -12a_0{a}_1{b}_0^2\lambda +15\alpha \mu a_1{a}_2{b}_0^2\lambda -52\mu a_1{a}_2{b}_0^2\lambda +3\alpha \mu a_0{a}_1{b}_1^2\lambda -12\mu a_0{a}_1{b}_1^2\lambda \hfill \\ \hfill & +20\mu a_0{b}_0{b}_1^2\lambda -2a_0{b}_0{b}_1^2\lambda -8\alpha a_0{a}_1^2{b}_0\lambda +12a_0{a}_1^2{b}_0\lambda -8\alpha a_0^2{a}_2{b}_0\lambda +12a_0^2{a}_2{b}_0\lambda \hfill \\ \hfill & -20\mu a_1{b}_0^2{b}_1\lambda +2a_1{b}_0^2{b}_1\lambda +8\alpha a_0^2{a}_1{b}_1\lambda -12a_0^2{a}_1{b}_1\lambda -3\alpha a_0^2{b}_0{b}_1\lambda +12a_0^2{b}_0{b}_1\lambda \hfill \\ \hfill & -3\alpha \mu a_1^2{b}_0{b}_1\lambda +12\mu a_1^2{b}_0{b}_1\lambda -6\alpha \mu a_0{a}_2{b}_0{b}_1\lambda +16\mu a_0{a}_2{b}_0{b}_1\lambda +8\mu a_1{b}_0^3+a_1{b}_0^3\hfill \\ \hfill & +\omega a_2{b}_0^3-8{\mu }^2{a}_0{b}_1^3-\mu a_0{b}_1^3+6\alpha a_0^2{a}_1^2-6a_0^2{a}_1^2+4\alpha \mu a_1^2{b}_0^2-14\mu a_1^2{b}_0^2-a_1^2{b}_0^2\hfill \\ \hfill & +6\alpha {\mu }^2{a}_2^2{b}_0^2-20{\mu }^2{a}_2^2{b}_0^2+8\alpha \mu a_0{a}_2{b}_0^2-32\mu a_0{a}_2{b}_0^2-2a_0{a}_2{b}_0^2+4\alpha \mu a_0^2{b}_1^2-14\mu a_0^2{b}_1^2\hfill \\ \hfill & -a_0^2{b}_1^2-2{\mu }^2{a}_0{a}_2{b}_1^2+3\omega a_0{b}_0{b}_1^2+8{\mu }^2{a}_1{b}_0{b}_1^2+\mu a_1{b}_0{b}_1^2+4\alpha a_0^3{a}_2-4a_0^3{a}_2-4\alpha \mu a_1^3{b}_0\hfill \\ \hfill & +6\mu a_1^3{b}_0-4\alpha a_0^2{a}_1{b}_0+6a_0^2{a}_1{b}_0-24\alpha \mu a_0{a}_1{a}_2{b}_0+36\mu a_0{a}_1{a}_2{b}_0+4\alpha a_0^3{b}_1-6a_0^3{b}_1\hfill \\ \hfill & +4\alpha \mu a_0{a}_1^2{b}_1-6\mu a_0{a}_1^2{b}_1-8\mu a_0{b}_0^2{b}_1-a_0{b}_0^2{b}_1+3\omega a_1{b}_0^2{b}_1+4{\mu }^2{a}_2{b}_0^2{b}_1\hfill \\ \hfill & +5\mu a_2{b}_0^2{b}_1+4\alpha \mu a_0^2{a}_2{b}_1-6\mu a_0^2{a}_2{b}_1-8\alpha \mu a_0{a}_1{b}_0{b}_1\hfill \\ \hfill & +28\mu a_0{a}_1{b}_0{b}_1-4a_0{a}_1{b}_0{b}_1+2{\mu }^2{a}_1{a}_2{b}_0{b}_1=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill e^{7Q\left(z\right)}:& a_1{b}_0^3{\lambda }^3-a_0{b}_0^2{b}_1{\lambda }^3+14\mu a_2{b}_0^3{\lambda }^2+\alpha a_0{a}_1{b}_0^2{\lambda }^2-4a_0{a}_1{b}_0^2{\lambda }^2+10\mu a_0{b}_0{b}_1^2{\lambda }^2\hfill \\ \hfill & -10\mu a_1{b}_0^2{b}_1{\lambda }^2-\alpha a_0^2{b}_0{b}_1{\lambda }^2+4a_0^2{b}_0{b}_1{\lambda }^2+8\mu a_1{b}_0^3\lambda +a_1{b}_0^3\lambda -12{\mu }^2{a}_0{b}_1^3\lambda \hfill \\ \hfill & +3\alpha \mu a_1^2{b}_0^2\lambda -10\mu a_1^2{b}_0^2\lambda +6\alpha \mu a_0{a}_2{b}_0^2\lambda -24\mu a_0{a}_2{b}_0^2\lambda +5\alpha \mu a_0^2{b}_1^2\lambda \hfill \\ \hfill & -18\mu a_0^2{b}_1^2\lambda +12{\mu }^2{a}_1{b}_0{b}_1^2\lambda -4\alpha a_0^2{a}_1{b}_0\lambda +6a_0^2{a}_1{b}_0\lambda +4\alpha a_0^3{b}_1\lambda -6a_0^3{b}_1\lambda \hfill \\ \hfill & -8\mu a_0{b}_0^2{b}_1\lambda -a_0{b}_0^2{b}_1\lambda -12{\mu }^2{a}_2{b}_0^2{b}_1\lambda -8\alpha \mu a_0{a}_1{b}_0{b}_1\lambda +28\mu a_0{a}_1{b}_0{b}_1\lambda \hfill \\ \hfill & +\omega a_1{b}_0^3+16{\mu }^2{a}_2{b}_0^3+2\mu a_2{b}_0^3+2\alpha \mu a_0{a}_1{b}_0^2-8\mu a_0{a}_1{b}_0^2-2a_0{a}_1{b}_0^2+6\alpha {\mu }^2{a}_1{a}_2{b}_0^2\hfill \\ \hfill & -20{\mu }^2{a}_1{a}_2{b}_0^2+2\alpha {\mu }^2{a}_0{a}_1{b}_1^2-8{\mu }^2{a}_0{a}_1{b}_1^2+8{\mu }^2{a}_0{b}_0{b}_1^2-2\mu a_0{b}_0{b}_1^2+4\alpha a_0^3{a}_1\hfill \\ \hfill & -4a_0^3{a}_1-8\alpha \mu a_0{a}_1^2{b}_0+12\mu a_0{a}_1^2{b}_0-8\alpha \mu a_0^2{a}_2{b}_0+12\mu a_0^2{a}_2{b}_0+3\omega a_0{b}_0^2{b}_1\hfill \\ \hfill & -8{\mu }^2{a}_1{b}_0^2{b}_1+2\mu a_1{b}_0^2{b}_1+8\alpha \mu a_0^2{a}_1{b}_1-12\mu a_0^2{a}_1{b}_1-2\alpha \mu a_0^2{b}_0{b}_1+8\mu a_0^2{b}_0{b}_1\hfill \\ \hfill & -2a_0^2{b}_0{b}_1-2\alpha {\mu }^2{a}_1^2{b}_0{b}_1+8{\mu }^2{a}_1^2{b}_0{b}_1-4\alpha {\mu }^2{a}_0{a}_2{b}_0{b}_1+12{\mu }^2{a}_0{a}_2{b}_0{b}_1=0,\hfill \end{array}$$ $$\begin{array}{cc}\hfill e^{8Q\left(z\right)}:& \alpha a_0^4-\alpha a_0^2{b}_0{b}_1\lambda \mu +\alpha a_1{a}_0{b}_0^2\lambda \mu +3\alpha a_0^2{b}_1^2{\mu }^2+2\alpha a_2{a}_0{b}_0^2{\mu }^2-4\alpha a_1{a}_0{b}_0{b}_1{\mu }^2\hfill \\ \hfill & +\alpha a_1^2{b}_0^2{\mu }^2+4\alpha a_0^3{b}_1\mu -4\alpha a_1{a}_0^2{b}_0\mu -a_0{b}_0^2{b}_1{\lambda }^2\mu +a_1{b}_0^3{\lambda }^2\mu +6a_0{b}_0{b}_1^2\lambda {\mu }^2\hfill \\ \hfill & +6a_2{b}_0^3\lambda {\mu }^2-6a_1{b}_0^2{b}_1\lambda {\mu }^2+4a_0^2{b}_0{b}_1\lambda \mu -4a_1{a}_0{b}_0^2\lambda \mu -6a_0{b}_1^3{\mu }^3+6a_1{b}_0{b}_1^2{\mu }^3\hfill \\ \hfill & -6a_2{b}_0^2{b}_1{\mu }^3-11a_0^2{b}_1^2{\mu }^2-8a_2{a}_0{b}_0^2{\mu }^2-2a_0{b}_0^2{b}_1{\mu }^2+14a_1{a}_0{b}_0{b}_1{\mu }^2+2a_1{b}_0^3{\mu }^2\hfill \\ \hfill & -3a_1^2{b}_0^2{\mu }^2-6a_0^3{b}_1\mu +6a_1{a}_0^2{b}_0\mu -a_0{b}_0^2{b}_1\mu +a_1{b}_0^3\mu +a_0{b}_0^3\omega -a_0^2{b}_0^2-a_0^4=0.\hfill \end{array}$$

References

1. Keller, E.F.; Segel, L.A. Traveling bands of chemotactic bacteria: A theoretical analysis. J. Theor. Biol.; 1971; 30, pp. 235-248. [DOI: https://dx.doi.org/10.1016/0022-5193(71)90051-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/4926702]

2. Arumugam, G.; Tyagi, J. Keller-Segel Chemotaxis Models: A Review. Acta Appl. Math.; 2020; 171, 6. [DOI: https://dx.doi.org/10.1007/s10440-020-00374-2]

3. Bhaya, D.; Takahashi, A.; Grossman, A.R. Light regulation of type IV pilus-dependent motility by chemosensor-like elements in Synechocystis PCC6803. Proc. Natl. Acad. Sci. USA; 2001; 98, pp. 7540-7545. [DOI: https://dx.doi.org/10.1073/pnas.131201098] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11404477]

4. Varuni, P.; Menon, S.N.; Menon, G.I. Phototaxis as a collective phenomenon in Cyanobacterial colonies. Sci. Rep.; 2017; 7, 17799. [DOI: https://dx.doi.org/10.1038/s41598-017-18160-w] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29259320]

5. Levy, D.; Requeijo, T. Modeling group dynamics of phototaxis: From particle systems to PDEs. Discret. Contin. Dyn. Syst.-B; 2008; 9, pp. 103-128. [DOI: https://dx.doi.org/10.3934/dcdsb.2008.9.103]

6. Levy, D.; Requeijo, T. Stochastic models for phototaxis. Bull. Math. Biol.; 2008; 70, pp. 1684-1706. [DOI: https://dx.doi.org/10.1007/s11538-008-9314-5]

7. Ha, S.; Levy, D. Particle, kinetic and fluid models for phototaxis. Discret. Contin. Dyn. Syst.-B; 2009; 12, pp. 77-108. [DOI: https://dx.doi.org/10.3934/dcdsb.2009.12.77]

8. Galante, A.; Wisen, S.; Bhaya, D.; Levy, D. Modeling local interactions during the motion of cyanobacteria. J. Theor. Biol.; 2012; 309, pp. 147-158. [DOI: https://dx.doi.org/10.1016/j.jtbi.2012.06.013]

9. Galante, A.; Levy, D. Modeling selective local interactions with memory. Phys. D Nonlinear Phenom.; 2013; 260, pp. 176-190. [DOI: https://dx.doi.org/10.1016/j.physd.2012.10.010] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24244060]

10. Weinberg, D.; Levy, D. Modeling selective local interactions with memory: Motion on a 2d lattice. Phys. D Nonlinear Phenom.; 2014; 278–279, pp. 13-30. [DOI: https://dx.doi.org/10.1016/j.physd.2014.04.001]

11. Drescher, K.; Goldstein, R.; Tuval, I. Fidelity of adaptive phototaxis. Proc. Natl. Acad. Sci. USA; 2010; 107, pp. 11171-11176. [DOI: https://dx.doi.org/10.1073/pnas.1000901107] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20534560]

12. Giometto, A.; Altermatt, F.; Maritan, A.; Stocker, R.; Rinaldo, A. Generalized receptor law governs phototaxis in the phytoplankton Euglena gracilis. Proc. Natl. Acad. Sci. USA; 2015; 112, pp. 7045-7050. [DOI: https://dx.doi.org/10.1073/pnas.1422922112]

13. Dervaux, J.; Resta, M.C.; Brunet, P. Light-controlled flows in active fluids. Nat. Phys.; 2017; 13, pp. 306-312. [DOI: https://dx.doi.org/10.1038/nphys3926]

14. Chavy-Waddy, P.; Kolokolnikov, T. A local PDE model of aggregation formation in bacterial colonies. Nonlinearity; 2016; 29, 3174. [DOI: https://dx.doi.org/10.1088/0951-7715/29/10/3174]

15. Bernoff, A.J.; Topaz, C.M. Biological aggregation driven by social and environmental factors: A nonlocal model and its degenerate Cahn–Hilliard approximation. SIAM J. Appl. Dyn. Syst.; 2016; 15, pp. 1528-1562. [DOI: https://dx.doi.org/10.1137/15M1031151]

16. Taranets, R.; Chugunova, M. Longtime dynamics of the PDE model for the motion toward light of bacterial colonies. Nonlinearity; 2018; 31, 887. [DOI: https://dx.doi.org/10.1088/1361-6544/aa999d]

17. Leptos, K.C.; Chioccioli, M.; Furlan, S.; Pesci, A.I.; Goldstein, R.E. Phototaxis of chlamydomonas arises from a tuned adaptive photoresponse shared with multicellular volvocine green algae. Phys. Rev. E; 2023; 107, 014404. [DOI: https://dx.doi.org/10.1103/PhysRevE.107.014404]

18. Ali, K.K.; Osman, M.; Baskonus, H.M.; Elazab, N.; Ilhan, E. Analytical and numerical study of the HIV-1 infection of CD4+ T-cells conformable fractional mathematical model that causes acquired immunodeficiency syndrome (AIDS) with the effect of antiviral drug therapy. Math. Methods Appl. Sci.; 2020; 46, pp. 7654-67670. [DOI: https://dx.doi.org/10.1002/mma.7022]

19. Kumar, D.; Seadawy, A.R.; Joardar, A.K. Modified Kudryashov method via new exact solutions for some conformable fractional differential equations arising in mathematical biology. Chin. J. Phys.; 2018; 56, pp. 75-85. [DOI: https://dx.doi.org/10.1016/j.cjph.2017.11.020]

20. Lee, D.; Huh, J.; Jeong, D.; Shin, J.; Yun, A.; Kim, J. Physical, mathematical, and numerical derivations of the Cahn–Hilliard equation. Comput. Mater. Sci.; 2014; 81, pp. 216-225. [DOI: https://dx.doi.org/10.1016/j.commatsci.2013.08.027]

21. Murray, J.D. Mathematical Biology I. An Introduction, Volume 17 of Interdisciplinary Applied Mathematics; Springer: Berlin/Heidelberg, Germany, 2002.

22. Murray, J.D. Mathematical Biology II: Spatial Models and Biomedical Applications, Volume 18 of Interdisciplinary Applied Mathematics; Springer: Berlin/Heidelberg, Germany, 2003.

23. Couzin, I.D.; Krause, J.; James, R.; Ruxton, G.D.; Franks, N.R. Collective memory and spatial sorting in animal groups. J. Theor. Biol.; 2002; 218, pp. 1-11. [DOI: https://dx.doi.org/10.1006/jtbi.2002.3065] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12297066]

24. Mogilner, A.; Edelstein-Keshet, L.; Bent, L.; Spiros, A. Mutual interactions, potentials, and individual distance in a social aggregation. J. Math. Biol.; 2003; 47, pp. 353-389. [DOI: https://dx.doi.org/10.1007/s00285-003-0209-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/14523578]

25. Volpert, V.A.; Volpert, A.I. Application of the Leray-Schauder method to the proof of the existence of wave solutions of parabolic systems. Sov. Math.; 1988; 37, pp. 138-141.

26. Volpert, A.I.; Volpert, V.A.; Volpert, V.A. Traveling Wave Solutions of Parabolic Systems. Translations of Mathematical Monographs; American Mathematical Society: Providence, RI, USA, 1994.

27. Kudryashov, N.A. One method for finding exact solutions of nonlinear differential equations. Commun. Nonlinear Sci. Numer. Simul.; 2012; 17, pp. 2248-2253. [DOI: https://dx.doi.org/10.1016/j.cnsns.2011.10.016]

28. Kaplan, M.; Bekir, A.; Akbulut, A. A generalized kudryashov method to some nonlinear evolution equations in mathematical physics. Nonlinear Dyn.; 2016; 85, pp. 2843-2850. [DOI: https://dx.doi.org/10.1007/s11071-016-2867-1]

29. Akbar, M.A.; Ali, N.H.M. Solitary wave solutions of the fourth order Boussinesq equation through the exp(-ϕ(η))-expansion method. SpringerPlus; 2014; 3, 344. [DOI: https://dx.doi.org/10.1186/2193-1801-3-344]

30. Uddin, S.; Alam, N.; Hossain, S.M.S.; Hasan, S. Some New Exact Traveling Wave Solutions to the (3+1)-Dimensional Zakharov-Kuznetsov Equation and the Burgers Equations via Exp(-ϕ(η))-Expansion Method. Front. Math. Its Appl.; 2014; 1, pp. 1-8.

31. Hafez, M.; Alam, M.N.; Akbar, M.A. Traveling wave solutions for some important coupled nonlinear physical models via the coupled Higgs equation and the Maccari system. J. King Saud Univ. Sci.; 2015; 27, pp. 105-112. [DOI: https://dx.doi.org/10.1016/j.jksus.2014.09.001]

32. He, J.; Wu, X. Exp-function method for nonlinear wave equations. Chaos Solitons Fractals; 2006; 30, pp. 700-708. [DOI: https://dx.doi.org/10.1016/j.chaos.2006.03.020]

33. Bulut, H. Application of the modified exponential function method to the Cahn-Allen equation. AIP Conf. Proc.; 2017; 1798, 020033.