Existence of the solitary wave solutions supported by the modiﬁed FitzHugh–Nagumo system

. We study a system of nonlinear differential equations simulating transport phenomena in active media. The model we are interested in is a generalization of the celebrated FitzHugh– Nagumo system describing the nerve impulse propagation in axon. The modeling system is shown to possesses soliton-like solutions under certain restrictions on the parameters. The results of theoretical studies are backed by the direct numerical simulation.


Introduction
Studies of traveling wave (TW) solutions to nonlinear evolution equations attract attention of many researchers. Such interest is quite natural due to the fact that the TW solutions play an important role in the description of nonlinear phenomena in various fields of natural sciences, such as combustion and detonation [25], mathematical biology [8,13,20,24], nonlinear optics [11], and hydrodynamics [5]. A significant achievement of the theory of nonlinear waves was the development of the theory of solitons at the second half of the XX century [5]. Stability and particle-like properties of solitons, exposed during interactions, are often attributed to the Hamiltonian nature and complete integrability of the corresponding equations manifested in presence of infinite hierarchy of conservation laws. However, there is a large number of evolutionary equations of the dissipative type, which also have soliton-type solutions. In contrast to the Hamiltonian systems, dissipative models possess soliton-like solutions only for the selected values of the parameters. Nevertheless, such solutions very often not only exhibit stability features, but also possess attracting properties [2,16] and are therefore of particular interest as potential carriers of energy and information in open dissipative systems.
The subject of this work is a study of the solitary wave solutions supported by the modification of the FitzHugh-Nagumo equations [8,20] taking into account the effects of relaxation. For the FitzHugh-Nagumo model, which is much simpler than the original Hodgkin-Huxley equations [13], it has become possible to show the existence of solutions describing moving pulses, analyze their stability, and confirm the existence of a threshold energy value below which the localized solutions are unstable. It should be noted that rigorous studies of the FitzHugh-Nagumo model have been and remain a challenge to date. Therefore, much simple model is put forward by McKean [15] in 1970. In subsequent years, it has been studied in detail [7,21]. Simplification is achieved in this model by replacing the cubic nonlinear function, present in the FitzHugh-Nagumo model, with a piecewise linear Z-shaped function. Later on, it has been proposed the following modification of the McKean system [17]: (1) w t = bv − dw.
A concept leading to the equation with τ > 0 is presented in paper [14]. The formal derivation of equation (1) is based on the replacement of the Fick's law stating the generalized thermodynamical flow-force relation in the balanced equation for the variable v with the Cattaneo's generalization [4] τ J t (t, x) + J(t, x) = −K∇Q(t, x), which takes into account the effects of memory connected with the presence of internal structure on mesoscale.
In the present work, we study the following system: where Substantiation for this type of models was first proposed in papers [6,18]. Let us note that the global existence and uniqueness results to a class of systems more general than (3) have been presented recently in papers [22,23]. Replacing a piecewise linear function with a function having the cubic nonlinearity greatly complicates the study since it eliminates the possibility of constructing exact solutions of the soliton type. There are two main trends of research of soliton solutions in FitzHugh-Nagumo-type systems. In the concept based on the so-called slow-fast systems' approach [3], the presence of a small parameter in the right side of the kinetic equation is essentially used. The approach we follow in this work is not directly related to the presence of a small parameter in the system, although the most important results are proved precisely for small values of the parameter , and so far there is no reason to say that they can be transferred to a more general case. In this approach [12], the phase trajectories of a dynamical system associated with the initial system of PDEs are considered, and their dependence on the parameters of the system are studied. The main problem we address in this work is the proof of existence of the solitary waves among the set of TW solutions of the modified FitzHugh-Nagumo system. The TW solutions in this case satisfy a system of ordinary differential equations, called the factorized system, which can be presented in the form of a multi-dimensional dynamical system. Solitary wave solutions are represented by the homoclinic loops, i.e. by the trajectories bi-asymptotic to a stationary point.
The structure of this work is following. In Section 2, we present a general model describing transport phenomena in active systems, manifesting effects of temporal nonlocality, and derive the conditions for consistency of the obtained model with system (3) that we are going to analyze. In Section 3, we prove the existence of the solitary wave solutions. The proof is based on a number of additional statements presented in the form of lemmas and propositions. The main result (Theorem 3) is stated at the end of Section 3. Along with this, the profiles of the solitary waves obtained numerically are presented. In Section 4, the results obtained are summarized, and the areas for further research are outlined.
2 A model of transport processes in active media, manifesting the effects of memory A hyperbolic modification of the FitzHugh-Nagumo equation can be obtained directly from balance equations in which memory effects are taken into account. Suppose that the function v(t, x) describes the density of a certain physical quantity. Then changes of this quantity inside a bounded set Ω ⊂ R n , n = 1, 2, 3, are given by the equation where q(t, ξ) is the flow of the considered physical quantity at the point ξ ∈ ∂Ω, n(ξ) is the unit vector perpendicular to the surface ∂Ω in the point ξ, directed outward, I[v(t, x), w(t, x)] is the source term, w(t, x) is the parameter satisfying the kinetic equation If the medium is close to equilibrium, the generalized thermodynamical flow q(t, ξ) and the generalized thermodynamical force −const · ∇v(t, ξ), where ∇ is the gradient operator, satisfy the equation similar to (2). In more general case, when memory effects associated, for example, with the presence of an internal structure are taken into account, the following relation is used [4,14]: where K(t, t ) is the kernel of relaxation. The most commonly used is the Maxwell-Cataneo kernel taking the form The parameter τ is called the time of relaxation. Inserting (6) with kernel (7) into equation (4), we get Differentiating equation (8) with respect to time variable and then replacing the surface integral in the r.h.s with the voluminal integral, we obtain where ∆ is the Laplace operator. From the arbitrainess of Ω it follows that this relation should be satisfied at an arbitrary point x. Thus, we get the following system: In paper by Engelbrecht [6], analogous system for one spatial coordinate is obtained from the phenomenological considerations. The author constructs a model governing the electric impulses transmission along the squid's axon. In contrast to the preceding works (cf., e.g., with [24]), he takes into account a self-inductance of axon, incorporating this way the effects of temporal nonlocality. Initial system from paper [6] is as follows: where ρ is the axon radius; R is the specific resistance; C a , L are, correspondingly, the specific capacity and inductivity; i a is the current; v is the voltage, I(v, w) is the source term, w is the internal variable satisfying the kinetic equation (5). We are going to exclude the variable i a from system (11)- (12). For this purpose, we differentiate equation (11) with respect to t and then multiply it by L/(πρ 2 ); next, we differentiate equation (12) with respect to x. As a result, we get the following system: Subtracting the second equation form the first one and adding equation (5), we obtain the system analogous to (9)-10: , we want to choose the function I(v, w) in such a way that the right-hand side of equation (13) coincide with that of equation (3). We assume that the function I(v, w) has the following form: Equating this function to f (v)−w, we determine the constants: Finally, introducing the scaling (1/ √ RC a ) ∂/∂x → ∂/∂x, we get the system Equation (14) differs from the first equation of system (3) by the presence of an additional factors at the term v t . As it was noted in the introduction, the main purpose of this paper is to rigorously prove the existence of a soliton solution for a hyperbolic analogue of the FitzHugh-Nagumo system. Due to the complexity of the problem, we restrict ourselves to proving the existence of such solutions for the model system (3). It is worth noting that numerical experiments confirm the presence of soliton regimes among the solutions of system (14)- (15), however, the rigorous proof of their existence is beyond the scope of this paper.

Existence of solitary wave solutions
In what follows, we are interested in the traveling wave solutions v(t, where c > 0 is the velocity of the traveling wave moving from right to left. Inserting these functions into (3), we get the system Making the substitutions β = 1/(1 − τ c 2 ), δ = /c and introducing new variable u = v , we obtain the following dynamical system: (we assume further on that β > 0). Our aim is to show that, on certain restrictions on the parameters' values, system (16) possesses homoclinic orbits corresponding to the soliton-like TW solutions of the initial system. We will consider a, τ , γ as auxiliary (fixed) parameters assuming that all of them are positive, whereas the parameters and c as the main ones. Note that, due to the assumption regarding the sign of the parameter β, the velocity c cannot be arbitrarily large (which is unphysical) because it should satisfy the inequality c 2 < 1/τ . We will also assume that The above restriction assures that (0, 0, 0) is the only stationary point of system (16).

Local invariant manifolds of the origin
Linearization matrix for system (16) at the origin takes the form The characteristic equation of the matrixÂ is as follows: Lemma 1. The matrixÂ has one positive eigenvalue λ 1 and a pair of eigenvalues λ 2,3 with negative real parts.
Remark 1. Let us note that the eigenvector corresponding to λ j takes the form Thus, under the restriction (17), there exists a one-dimensional local invariant unstable manifold W u loc tangent to the eigenvector Y 1 at the origin and a two-dimensional local invariant stable manifold W s loc tangent to the plane spanned by the vectors Y 2 , Y 3 . Note that W u loc consists of two branches, one is tangent to Y 1 and pointed into the first octant (which is of interest to us), while the other is tangent to −Y 1 .
3.2 Behavior of saddle separatrices in cases = 0 and 0 < c 1 In this subsection, we use the condition = 0, allowing to separate the first two equations from the third one, having the trivial solution w = const. Arguments that will be given later on allow us without loss of generality to restrict consideration to the case w = 0.
The remaining system then takes the form For c = 0, system (18) can be presented in the Hamiltonian form with the Hamiltonian This Hamiltonian corresponds to the so-called system with one degree of freedom. In most cases, such a system can be fully analyzed by qualitative methods (see e.g. [1]). System (18) has the following stationary points: A = (0, 0), B = (a, 0), and C = (1, 0). The first coordinate of each of the points corresponds to the extremal value of the potential energy The local minimum of (19) corresponds to the center, while the local maximum corresponds to the saddle point [1]. Analyzing the signs of the second derivatives of the function U p (v) at the corresponding points, one can conclude that for 0 < a < 1, the points A and C are saddles, while B is a center. When a > 1, the point A and B, in turn, are saddles, while the point C is a center. We are looking for conditions assuring the existence of the trajectory doubly asymptotic to the saddle point A. Such trajectory exists when either a < 1/2 and B is the center or when a > 2 and the center is located in the point C. It turns out that the case a > 2 is not independent. Indeed, using in the case a > 2 the scaling transformation T = a 2 t, X = ax, τ =τ /a 2 ,ã = a −1 < 1/2, v = v/a,w = w/a,˜ = /a 2 , one can write down the source system as follows: So we will assume from now on that 0 < a < 1/2. Under the given assumption, the stationary point B = (a, 0) is a center, and thus, there is an open set U containing B, In what follows, we will need the information about the behavior of the separatrices of the saddle point A. We assume that 0 < c 1 and denote the stable and unstable separatrices of the saddle point A located in the right half-plane by q s c (ξ) and q u c (ξ), correspondingly. The following statement holds true. Proposition 1. The saddle separatrice q u c (ξ) directed towards the first quadrant • intersects the horizontal axis at some point v * such that a < v * < 1; • q u c (ξ) = (v(ξ), u(ξ)) tends to (−∞, −∞) as ξ → ω, where 0 < ω +∞. Proof. The proof of the first item is based on the Melnikov theory [10,19]. Assuming that c is small, we can present system (18), up to O(c 2 ), in the following form: where F = (∂H/∂u, −∂H/∂v) = (u, −f (v)), G = (0, u). For c > 0, the stationary point B turns into unstable focus, and stable and unstable separatrices of the saddle A do not form a closed loop any more. We want to trace what happens with the stable and unstable separatrices when 0 < c is small. Let us denote the point at which the homoclinic loop corresponding to c = 0 intersects the horizontal axis by p. For 0 < c 1, the stable and unstable separatrices will be located in the neighborhood of the homoclinic curve. We denote by q s c (0) and q u c (0) the points at which the stable and unstable separatrices intersect the horizontal axis (see Fig. 1). Up to O(c 2 ), the projection of the vector q u c (0) − q s c (0) onto the vector F ⊥ = (f (v), u) is given by the Melnikov integral [10,19] Thus the stable and unstable separatrices in the case c > 0 form the configuration shown on Fig. 1. Further behavior of the saddle separatrice q u c (ξ) is following. After the intersection of the horizontal axis, it enters the fourth quadrant, and its coordinate u(ξ) remains negative further on since the trajectory is separated from the upper half-plane by the separatrices of the saddle points A and C. Thus the coordinate v(ξ) decreases as ξ grows. The coordinate u(ξ), in turn, decreases when v(ξ) > a and increases when 0 < v(ξ) < a remaining negative. It becomes decreasing function again when v(ξ) is negative, and from this instant q u c (ξ) monotonically tends to (−∞, −∞).
Let us analyze the behavior of solutions for nonzero and c satisfying the conditions 0 < c 1. Up to the term of the order O(c 2 ), system (16) in this case can be presented as follows: where δ = /c 1. We are interested in the behavior of the trajectory q u c (ξ) = (v u c (ξ), u u c (ξ), w u c (ξ)) being the three-dimensional deformation of the trajectory q u c (ξ) and satisfying the condition lim ξ→−∞ q u c (ξ) = 0. The analysis of the linearization of system (20) shows that such deformation does exist. Without the loss of generality, we can assume that v u Presenting the solution to the third equation of system (20) in the form w = δw 1 + O(δ 2 ), we obtain with the specified accuracy the following representation: where v(z) is the first coordinate of the unstable saddle separatrice of the stationary point A of system (18). For δ 1, v u c (ξ) still dominates the behavior of the third variable, and therefore w u c in the r.h.s of the second equation does not influence the qualitative behavior neither the variable u u c (ξ) nor the variable v u c (ξ), which becomes negative and monotonically decreasing from some instant. And when the function |v u c (ξ)| becomes large enough, all three functions monotonically tend to −∞. Let us formulate the result obtained as follows.

Corollary 1.
There exist an open set ∆ in the space of the parameters (c, ), placed at the first quadrant and adjacent to the horizontal axis, such that for all (c, ) ∈ ∆, the phase trajectory q u c (ξ) satisfies the condition lim ξ→+∞ q u c (ξ) = (−∞, −∞, −∞).

The sets positively invariant with respect to the phase flow of the factorized system
In this subsection, we will prove the following lemma. Proof. In the proof below, as well as in the proofs of the subsequent assertions, we mainly follow the plan drawn in [12]. Thus, suppose that E + is not positively invariant with respect to the φ t , and the solution q(ξ) = (v(ξ), u(ξ), w(ξ)) satisfying q(0) ∈ E + leaves the set E + for the first time at ξ 1 > 0, so one of the features characterizing this set fails. To begin with, let us observe that the equality v(ξ 1 ) = 1 cannot be true since v(·) is growing on the segment (0, ξ 1 ) and v(0) > 1. For the same reason, u(ξ 1 ) cannot be equal to zero. Now let us address the function w (·). The relation w | (0,ξ1) > 0, together with the supposition w (ξ 1 ) = 0, imply the inequality w (ξ 1 ) 0, but w (ξ 1 ) = ( /c)u(ξ 1 ) > 0, hence we get the contradiction. Now, let us consider u (ξ). In accordance with the above assumptions, u (0) = β{cu(0) + w(0) − f [v(0)]} > 0, and since u(·), w(·), and v(·) are growing functions on the segment (0, ξ 1 ), and so is −f [v(·)], then u (ξ 1 ) cannot be zero as well. The positive invariance of the set E − is shown just in the same manner.

Asymptotic behavior of the unstable invariant manifold
We still assume that q u c, (ξ) = (v u c, (ξ), u u c, (ξ), w u c, (ξ)) is the unstable invariant manifold of the stationary point (0, 0, 0) corresponding to the given values of the parameters c, . However, for simplicity, from now on we will omit the superscript. Without the loss of generality, we will also assume that v c, (0) = a and u c, (ξ) > 0 when ξ 0. Let us define the following subsets of the set Ω = {(c, ): 0 < c < 1/ √ τ , 0} : ∈ Ω: q c, (ξ) is bounded , where 0 < ω i +∞, i = 0, 1. We will show that there is no other possible behavior of the trajectory q c, (ξ) differing from that presented above.
Lemma 3. The following statement is true: Proof. We construct such a rectangle A in the plane (v, w) that there will be only three possibilities: • Solution is bounded by this rectangle; • Solution leaves the rectangle and after that it enters E + , which implies that q c, (ξ) → (+∞, +∞, +∞); • Solution leaves the rectangle and after that it enters E − , which implies that q c, (ξ) → (−∞, −∞, −∞).
Just in the similar way, we can prove that if w > v/γ, then the phase trajectory approaches the left boundary of the rectangle A and ultimately falls into the set E − , so the corresponding pair (c, ) belongs to the set Ω 3 .
To complete the proof of this part, it suffices to note that the phase trajectory cannot leave the rectangle by its top or bottom. But it is quite evident since, at the points belonging to the upper and lower boundaries, the vector field is directed inward the rectangle A. And if the solution q c, does not leave the set A, then the pair (c, ) belongs to Ω 1 , and since there is no other choice, the statement is completely proved.

Insight into the structure of the subsets Ω 2
In this subsection, the geometry of subset of the parameters (c, ) ∈ Ω for which all the components of the vector-function q u c, (ξ) go to +∞ will be highlighted. Theorem 1. There exists values c 1 > 0 and 1 > 0 such that: Proof. Without the loss of generality, we assume that v(0) = a > 0, u(0) > 0, and w(0) > 0. The positivity of all components of the vector q(0) appears from the fact that for ξ −1, the phase trajectory is close to the eigenvector Y 1 = (1, λ 1 , /(c λ 1 + γ)), and for ξ ∈ (−∞, 0), all three components are nondecreasing. Indeed, for ξ −1, w(ξ)/v(ξ) ≈ /(cλ 1 + γ), and if the r.h.s is less that γ −1 , then the projection of q(ξ) onto the plane (v, w) lies below the line w = v/γ. But this requirement is equivalent to the inequality γ < cλ 1 + γ, which is true because all the parameters are positive. For ξ ∈ (−∞, 0), the r.h.s of the third equation of system (16) remains nonnegative because, at the instant when the projection of the trajectory onto the plane (v, w) approaches the line w = v/γ, v(ξ) is positive, and the projection cannot cross the line w = v/γ. This, in turn, implies that all components of the vector q(ξ) for ξ ∈ (−∞, 0) are positive, and besides v(0) − γw(0) 0. Further, as long as v(ξ) is an increasing function, we can assume that u(t) = U (v(t)), w(t) = W (v(t)). Then The r.h.s of equation (23) is nonnegative until U is positive. But U (a) > 0, and under the above supposition, is nondecreasing as v > a, and q(ξ) attains the set E + . Now let us assume that 0 < U (a) f max /c and W (a) > f max . Then we get the estimation and q(ξ) attains the set E + . Now it is necessary to find the conditions assuring that the inequality W (a) > f max is fulfilled. Under the assumption U (a) f max /c, the inequality takes place on the segment v ∈ (0, a). Applying the substitution W (v) = A(v)e −ρv , ρ = γ/f max , we get the inequality which, after the integration w.r.t. v on the segment (0, a), takes the form where θ = /f max . From this we get the inequality where 0 < µ < 1. So, if > 1 = 2f max /(µa 2 ), then, regardless of the value of c > 0, W (a) > f max . Now let us estimate c 1 . Equation (22) can be rewritten in the form from which appears the inequality .
To complete the proof, we just show that the inequalities c > f max /H(a) and c < 1/ √ τ are compatible. The first one is equivalent to so the inequalities are compatible, and the statement is completely proved.
Now, let us show that the following assertion is true.
located between two parabolas, then the only bounded solution possible is the trivial solution q c, (ξ) = 0.
Proof. Suppose that, under the above conditions, there exists a nontrivial bounded solution q c, (ξ). By analogy with [12] we consider the function which, under certain conditions, is monotonically decreasing on the solutions of system (16). Differentiating the function G(v, u, w) and taking into account (16), we obtain nullifies the coefficients at the nondiagonal quadratic terms. Assuming that E = −2cβ, we get the following expression: In view of assumption (17), expression (26) will be negative if the parameters c, belong to the set Φ and v 2 +w 2 = 0. So G[v(ξ), u(ξ), w(ξ)] should tend monotonically to a finite value differing from zero as ξ tends to +∞, but this is impossible since the origin is the only stationary point of system (16).  Proof. All the solutions q c, (ξ) corresponding to (c, ) ∈ Φ are unbounded. On the other hand, the set Φ is open and has nonempty intersection with the set Ω 2 (see Fig. 3(a)).
It is seen from the geometry of the open sets Ω 2 and the set ∆ ∈ Ω 3 , adjacent to the horizontal axis on the basis of the Corollary 1 (see Fig. 3(a)), that there should exist two subsets of the set Ω 1 lying between them. One of these subsets is adjacent to the origin, while the second one lies closer to the line c = c 1 . So the line = const for sufficiently small > 0 intersects the set Ω 1 at least in two points separated from each other. It remains to show that there are pairs (c, ) corresponding to the orbits bi-asymptotic to the origin.

Solutions corresponding to the solitary waves
Let us consider the set where v min is the point of a local minimum of the function f (v) on the segment (0, 1), and Λ = (c, ) ∈ Ω: solution q c, (ξ) intersect Σ exactly two times and after that does not intersects the region v v min .
It follows from the definition of the set Λ that Ω 2 ∩Λ = ∅. On the other hand, the arguments following the proof of Proposition 1 imply that the set Λ contains a nonempty open subset of points belonging to Ω 3 . And now we are going to prove the following assertion.
Proof. The proof of this theorem is based on the ideas underlying the proof of the Lemma 9 of the paper [12], so we try to adhere to the notation that was adopted in this work.
Let us suppose the opposite, namely, that to every pair (c, ) belonging to the set Ω 3 ∩ ∂Λ corresponds the orbit q c, (ξ) with the following properties: • The orbit q c, (ξ) starts from the origin and points to the first octant as ξ −1. We can assume, without the loss of generality, that v (ξ) and u (ξ) are positive on the interval (−∞, 0) and v(0) = a; • At some value of the argument, say ξ = s − > 0, the orbit crosses the set Σ, for the first time, intersecting it at a point belonging to the plane v > v min , u = 0, and next, at ξ =š > s − , intersects Σ for the second time at a point belonging to the plane v = v min , u < 0; • Before crossing the plane v = 0 and going to infinity (suppose that such intersection take place at ξ = s 1 ), the orbit must touch the set Σ, say, at ξ = s 0 >š (otherwise it does not belong to the set ∂Λ). Analysis of the first equation of system (16) tells us that the touch point must be located at the intersection of planes {v = v min , u < 0} and {v > v min , u = 0}.
Looking at the second equation of system (16), we easily conclude that It will be shown below that, for sufficiently small , there does not exist the trajectory q c, with (c, ) ∈ Ω 3 ∩ ∂Λ, which is characterized by the relations v(s 0 ) = v min , u(s 0 ) = 0, u (s 0 ) 0.
where u = 0. From this appears that u (σ −1 ) 0. Evaluating the behavior of u , we will show that for sufficiently small , this is not true. Since u(σ 1 ) α, then we get at ξ = σ −1 the estimation Since the trajectory leaves the origin pointing towards the first octant, its projection onto the plane (v, w) should still be in the rectangle A shown in Fig. 2 when ξ = σ −1 (otherwise the phase trajectory cannot touch the set Σ any more). For this rectangle, the following estimation holds: |v−γw| 2Ṽ , whereṼ = max{|v 1 |, v 2 }. For the parameter δ = /c, we get the estimation , then we get the contradiction. And this proves the statement.  where F = (u, −f (v)), K = (0, cu + δw 1 ), and w 1 , up to O(δ 2 ), is given by formula (21). So the Melnikov integral [10] will take the form It is quite obvious that, at a certain ratio between c and , the right-hand side will be equal to zero. There is nothing surprising that the Melnikov method catches only one homoclinic trajectory for a fixed value of . This is due to the fact that for the value of the parameter c lying close to the line c = c 1 , the assertion 0 < c 1 may be false.
The conclusions of theoretical studies indicating the presence of a pair of homoclinic trajectories under the above restrictions on the parameters are verified using numerical simulation. Numerical experiments conducted at a = 0.13, = 0.006991097526935545, γ = 2.706215020212898, and τ = 14.554975027077534 confirm the existence of a pair of homoclinic solutions corresponding to c = 0.195747986 (slow solitary wave) and c = 0.225305407 (fast solitary wave). The corresponding points are located in the plane Ω to the right of the line = βN c 2 (see Fig. 4), which fully agrees with the results of theoretical analysis. Figure 4(b) shows the graphs of slow and fast solitons on the physical plane.

Conclusion
It was proved in this paper that system (3), under certain restrictions on the parameters, possesses a pair of soliton-like traveling wave solutions. These solutions correspond to the homoclinic trajectories of the factorized system (16). The existence of the solitary wave solutions is confirmed by numerical experiments, the results of which completely agree with the conclusions of the theoretical analysis of system (16). Our further efforts will be purposed at studying the stability of the solitary wave solutions found, their dynamic properties, as well as identifying the presence of multi-hump wave patterns and investigating their properties. Let us note that preliminary results concerning these issues can be seen in [9].