Homoclinic Chaos Unveiling Quorum Sensing Dynamics
Mariana Harris Pablo Aguirre Víctor F. Breña-Medina

A generalization of K-theory to operator systems
    Walter D. van Suijlekom
    September 4, 2024
================================================









§ ABSTRACT 
Quorum sensing orchestrates bacterial communication, which is vital for bacteria's population behaviour. We propose a mathematical model that unveils chaotic dynamics within quorum sensing networks, challenging predictability. 
The model considers the interaction between autoinducers (molecular signalling) and two subtypes of bacteria. We analyze the different dynamical scenarios to find parameter regimes for long-term steady-state behaviour, periodic oscillations, and even chaos. In the latter case, we find that the complicated dynamics can be explained by the presence of homoclinic Shilnikov bifurcations.
  


Quorum sensing modelling, Synchronisation, Shilnikov homoclinic bifurcation, chaos. 



34C23, 37G05, 37G15, 37C29, 37G20, 92B25.


myheadings
plain
M. Harris, P. Aguirre and V. Breña–Medina
Homoclinic Chaos and Quorum Sensing 







§ INTRODUCTION





Quorum sensing (QS) —also known as auto-induction— is a mechanism of regulation of gene expression which allows communication in a cell population <cit.>. In the case of bacteria, this process is a density dependent behaviour, which consists of cell-to-cell communication through the release of chemical signal molecules known as autoinducers. This communication enables a population to explore a given medium and alter its behaviour in response to its own fluctuation, which occurs from changes in the number of species present in the community. Specifically, the basis of the QS mechanism is the production and release of autoinducers, whose concentration increases as a function of bacterial population density. When the population density exceeds a threshold, the bacteria detect the autoinducers which then activate genes that switch on a behavioural trait. Autoinducers located in the external medium diffuse and bind with a specific protein within the bacteria. A protein-autoinducer complex is then formed that attaches to a region of the bacteria's DNA. This in turn regulates the production of the autoinducers enhancing a specific regulated cell-density behaviour; see, for instance, <cit.>.

The QS phenomenon was first observed by Nealson and  Hastings in 1970 in a bioluminescent bacteria called Vibrio fischeri <cit.>. This bacteria forms a mutually beneficial symbiotic relationship with some species of squid such as the Hawaiian squid Euprymna scolopes. The V. fischeri live in the squid's light organ which provides nourishment, allowing the bacteria to proliferate. When the bacteria reach a high population density the genes involved in bioluminescence are expressed. The produced light allows the squid to mask its shadow and avoid predation <cit.>. Following the discovery of the V. fisheri's density-dependent bioluminescent behaviour, other species of bacteria have been discovered to employ a QS mechanism; for instance, the bacteria Chromobacterium violeceum gives rise to the production of a purple pigment <cit.>. Also, biofilms may present challenges  and opportunities in the industry setting as these may cause blockages in specialised machinery <cit.>, gene expression <cit.>, disrupting QS signalling (a process known as quorum quenching) <cit.>, and inhibiting pathogen virulence <cit.>. 

The QS phenomena has been studied from either deterministic as well as stochastic and hybrid mathematical modelling approaches, which are inspired by many biologically plausible perspectives. For instance, in <cit.> a plant-pathogen population controlling virulence on leaf surfaces is investigated by means of a continuous-time Markov process, where a linear birth and a logistic-death-migration population processes along with an autocatalytic mechanism for acyl homoserine lactone autoinducers concentration are taken into account. Their findings include that QS mechanism and autoinducers diffusion follow an inversely proportional relationship. Another interesting work can be seen in <cit.>; there, extracellular polymeric substances (EPS) production is modelled in a growing biofilm under various environmental conditions, which yields to a reaction-diffusion system where the diffusion coefficient is density-dependent. The authors found that QS-induced EPS production permits a biofilm to switch from a colonisation state to a protection state, which is a crucial characteristic of QS. Even further, in <cit.> an anomalous diffusion process is taken into account, which captures a delayed maundering of substances as a consequence of transcriptional features. On the other hand, a deterministic gene regulatory network point of view is addressed in <cit.> by considering a QS pathway involving multiple feedforward and negative feedback loops aside from transcription time delays in a cancer drug realising scheme. A key result from their manuscripts lie on the appearance of Hopf bifurcations, other secondary oscillatory-induced bifurcations and a time-delay threshold, which coordinate self-sustained oscillating features. This small sample of approaches pays special attention to transcriptional time-delay as well as on-and-off activating crucial switches. Our area of focus is primarily directed towards the later. 

In this work we propose a simple model which consists of features of a bacterial population and its inherent autoinducer concentration interaction. In other words, we assume that the population of bacteria is locally activated and inhibited by the production of autoinducers.
In order to capture the key qualitative ingredients of QS mechanisms in bacteria, we follow a simplified approach. In so doing, we assume that the autoinducers bind to receptors that enhance the expression of a particular gene, and its production also directly depends on the bacteria population (e. g. <cit.>). When a bacterium responds to an external stimulus by increasing the amount of a cellular component, the process is up-regulated, and it is down-regulated otherwise. To put it another way, the dynamics of the bacteria are regulated by both, autoinducers as well as growth of the bacteria population itself. Hence, the population is considered to be composed of cell sub-populations in two different states: up-regulated and down-regulated bacteria, namely motile m and static s, respectively, and the concentration of autoinducers is given by q. As autoinducers are produced by both sub-populations at a rate r (see <cit.>), and the bacteria can either switch or remain in their category, we assume that:
[label=(*)]

  * the growth rate of the motile bacteria population follows an autocatalytic process and is proportional to the probability that motile bacteria switch on-and-off in their current category. The response rate of motile bacteria to autoinducer concentration is then given by k_1 and is up-regulated by its own production and down-regulated by the production of static bacteria; that is, the dynamics follows an activator-inhibitor-like behaviour, where m acts as the activator and s  as the inhibitor. We also assume that the media may be overpopulated, so γ is a saturation parameter;

  * the production of static bacteria is assumed as a result of a cross-catalytic reaction of motile bacteria; this interaction may be taken as an indirect down-regulation process, which is modulated by k_2;

  * the entire bacteria population is constantly produced at a rate (1+ε)α, where 0< ε <1 is the ratio between both bacteria sub-populations growth rates. The autoinducers, motile and static bacteria decay at rates μ_1, μ_2 and μ_3, respectively.  

Putting everything together, we obtain that the dynamical local interaction is governed by the system 


    X:{[                         q̇                          =              r(m+s)-μ_1q ,;                          ṁ                          = k_1qm^2s(1+γ m^2)+α-μ_2m ,;                          ṡ                          =          k_2qm^2+ϵα-μ_3s . ].


Essentially, model (<ref>) captures the dynamical behaviour of QS, namely, the concentration of autoinducers depends on the production by bacteria in both states. When this concentration lies beyond a certain threshold, bacteria are expected to react in a synchronised-like way, as has been reported in <cit.>, for instance. We have identified this feature by means of slowly varying a crucial parameter as can be seen further in section <ref>, and analyse the different dynamical scenarios of (<ref>) to find parameter regimes for steady state long term behaviour and the observed periodic oscillations. However, we also find chaos which can be explained by the presence of Shilnikov homoclinic bifurcations.
 
The present manuscript is organised as follows: throughout section <ref> we prove that the model solutions are always non-negative for feasible initial conditions. In section <ref> we take into consideration a constant autoinducer concentration and study the local stability of a reduced system. This is followed in  Section <ref> by a bifurcation analysis where the role of a key parameter triggers leading nonlinear events, namely Hopf and Bogdanov–Takens bifurcations. Normal forms for both dynamical phenomena are also computed. From there, a thorough numerical bifurcation analysis is performed in section <ref> in the full system, where conditions to the Shilnikov homoclinic chaos mechanism are found and numerically explored. Concluding remarks can be found in section <ref>.   




§ NON-NEGATIVE SOLUTIONS FOR REALISTIC INITIAL CONDITIONS




Model (<ref>) is well-posed in the sense that every solution starting from a realistic (positive) initial condition remains non-negative.  In what follows we see (<ref>) as a vector field defined in the set Ω:={(q,m,s)∈ℝ^3|  q≥0, m≥0, s>0} and use the following standard notation <cit.> for its components: X(q,m,s)=X_1(q,m,s)∂∂ q+X_2(q,m,s)∂∂ m+X_3(q,m,s)∂∂ s.

Let us study the behaviour of  X in ∂Ω to show that Ω is invariant.
Let us subdivide ∂Ω into the following coordinate planes: ∂Ω=ω_qm∪ω_qs∪ω_ms, where ω_qm={(q,m,s)∈ℝ^3|  q≥0, m≥0, s=0}, ω_qs={(q,m,s)∈ℝ^3|  q≥0, m=0, s>0}, and ω_ms={(q,m,s)∈ℝ^3|  q=0, m≥0, s>0}.

The restriction of (<ref>) to the plane ω_qs is X(q,0,s)=(rs-μ_1 q)∂∂ q+α∂∂ m-(ϵα+μ_3 s)∂∂ s. Since α>0, the vector field (<ref>) on ω_qs points towards the interior of Ω. Similarly, the restriction of (<ref>) to the plane ω_ms is X(0,m,s)=r(m+s)∂∂ q+(α-μ_2 m)∂∂ m-(ϵα+μ_3 s)∂∂ s. Since r>0, it follows that the vector field X in ω_ms points towards the interior of Ω.

On the other hand, system (<ref>) is not defined in the plane ω_qm. In order to analyse the behaviour of (<ref>) near ω_qm, let us consider the time scaling t↦ st. In this way, (<ref>) becomes


    X_0:{[                        q̇                         =            r(m+s)s-μ_1qs,;                         ṁ                         = k_1qm^21+γ m^2+α s-μ_2ms,;                         ṡ                         =     k_2qm^2s+ϵα s-μ_3s^2. ].


System (<ref>) is topologically equivalent to (<ref>) in Ω. Moreover, (<ref>) is well defined for s=0 and, hence, it can be continually extended to the boundary plane ω_qm⊂∂Ω. The restriction of (<ref>) to the plane ω_qm is X_0(q,m,0)=k_1qm^21+γ m^2∂∂ m. It follows that the (unique) solution of (<ref>) with initial condition (q(0),m(0),s(0))=(q_0,m_0,0)∈ω_qm is given by


    {[                                                              q(t)                                                                 =                                                              q_0,;                                                              m(t)                                                                 = -1+k_1q_0m_0t+γ m_0^2+√(4γ m_0^2+(k_1tq_0m_0+γ m_0^2-1)^2)2γ m_0,;                                                              s(t)                                                                 =                                                                0. ].


As a consequence, the set ω_qm is an invariant plane of (<ref>) which consists of a continuum of straight lines 
parallel to the m-axis parameterized by (<ref>). Moreover, both axes q=0 and m=0 consist of a continuum of equilibria.

Let us now study the behaviour of orbits of (<ref>) near the invariant plane ω_qm. Let us search for the solution of (<ref>) with initial condition (q(0),m(0),s(0))=(q_0,m_0,δ)∈ int(Ω), with 0<δ≪1 sufficiently small. Specifically, consider a solution of the form:


    {[                    q(t)                       = q_0(t)+δ q_1(t)+O(δ^2),;                    m(t)                       = m_0(t)+δ m_1(t)+O(δ^2),;                    s(t)                       = s_0(t)+δ s_1(t)+O(δ^2), ].


where (q_0(t),m_0(t),s_0(t)) is the solution of (<ref>) in the limit as δ→0 and, hence, it is given by (<ref>). In this way, (<ref>) is expressed as


    {[                    q(t)                       =    q_0+δ q_1(t)+O(δ^2),;                    m(t)                       = m_0(t)+δ m_1(t)+O(δ^2),;                    s(t)                       =        δ s_1(t)+O(δ^2), ].


with m_0(t)= -1+k_1q_0m_0t+γ m_0^2+√(4γ m_0^2+(k_1tq_0m_0+γ m_0^2-1)^2)2γ m_0. It follows that the higher order terms of (<ref>) must satisfy the initial conditions (q_1(0),m_1(0),s_1(0))=(0,0,1) for n=1, and (q_n(0),m_n(0),s_n(0))=(0,0,0) for n≥2. 
Substitution of (<ref>) into (<ref>) leads to the following differential equations for the O(δ) terms of q(t):


    δq̇_1(t)    =    r(m_0(t)+δ m_1(t)+…+δ s_1(t)+…)(δ s_1(t)+…)-μ_1(q_0+δ q_1(t)+…)(δ s_1(t)+…), 
        =    δ rm_0(t)s_1(t)-δμ_1 q_0s_1(t)+O(δ^2).


Hence,



    q̇_1(t)=rm_0(t)s_1(t)-δμ_1 q_0s_1(t),  q_1(0)=0.



Similarly,

    ṁ(t)=ṁ_0(t)+δṁ_1(t) +O(δ^2) = k_1q(t)m^2(t)1+γ m^2(t)+α s(t)-μ_2 m(t)s(t).

For 0<δ≪1 sufficiently small we have:


    ṁ_0(t)+δṁ_1(t) +O(δ^2)    =   k_1q_0m_0^2(t)1+γ m_0^2(t)+δ( k_1q_1(t)m_0^2(t)1+γ m_0^2(t)- 2γ k_1q_0m_0^3(t)m_1(t)(1+γ m_0^2(t))^2+ 2k_1q_0m_0^2(t)m_1(t)1+γ m_0^2(t))
     
           + δα s_1(t)-δμ_2 m_0(t)s_1(t)+O(δ^2).


Thus,


    ṁ_1(t)=k_1q_1(t)m_0^2(t)1+γ m_0^2(t)- 2γ k_1q_0m_0^3(t)m_1(t)(1+γ m_0^2(t))^2+ 2k_1q_0m_0^2(t)m_1(t)1+γ m_0^2(t)+ α s_1(t)-μ_2 m_0(t)s_1(t),   m_1(0)=0.


Likewise, 


    δṡ_1(t)    =    k_2(q_0+δ q_1(t)+…)(m_0(t)+δ m_1(t)+…)^2(δ s_1(t)+…)+ϵα(δ s_1(t)+…)
           - μ_3(δ s_1(t)+…)^2 = δ k_2 q_0m_0^2(t)s_1(t)+δϵα s_1(t)+O(δ^2).


Hence,


    ṡ_1(t)=(k_2 q_0m_0^2(t)+ϵα)s_1(t),  s_1(0)=1.



The initial value problem (<ref>) defines any solution of (<ref>) starting at a distance 0<δ≪1 from the invariant plane ω_qm with accuracy  O(1/δ). In particular, 
since k_2 q_0m_0^2(t)+ϵα>0, it follows from (<ref>) that s_1(t) is an increasing function for every t>0. Hence, the component s(t) in (<ref>) is increasing near the plane s=0. Therefore, no trajectory of (<ref>) with initial condition in int(Ω) can reach the boundary plane ω_qm in (finite or infinite) positive time.

Since (<ref>) is topologically equivalent to (<ref>) in Ω, we conclude from the analysis in this section that every trajectory of the original system (<ref>) starting in Ω remains non-negative.





§ CONSTANT AUTOINDUCER CONCENTRATION





To shed light on the understanding of QS we now analyse the bacteria interaction dynamics under the assumption of a constant autoinducer concentration q_0 > 0. In so doing, we have η_1 = k_1q_0 and η_2 = k_2q_0 as the motile and static bacteria response rates respectively. From (<ref>), upon substituting the rescaled variables for the bacteria population u= (η_2/η_1)m and v= (μ_3η_2/η_1^2)s as well as τ = μ_3 t for time, and the new parameters

    K = γ(η_1η_2)^2 ,     b = μ_2μ_3 ,     a= η_2 αη_1 μ_3 ,     e = μ_3 εη_1,

we obtain the following 2D constant autoinducer system  


    {[               u̇ =f(u,v), f(u,v)=u^2v(1+Ku^2)+a-bu,;              v̇ = g(u,v),          g(u,v)=u^2+ae-v. ].



Notice that key parameters arise so that: 
[label=(*)]

  * K holds a saturation role, 

  * b characterises the decaying bacteria rate, and 

  * a and e capture bacteria production-related roles as well as inversely and directly proportional dependence of the constant autoinducers' concentration and decaying rate variations, respectively. The latter shows that parameter product ae only depends on the bacteria response to autoinducers capacities and production of static subpopulation.









 §.§ Number of steady-states
 




As the understanding of interactions between motile and static bacteria is the backbone of communication mechanism via autoinducers, we seek the circumstances affecting the way subpopulations coexist or go extinct. Hence, we first look for positive steady-states. Consider that nullclines of (<ref>) satisfy


    v = u^2(1+Ku^2)(bu-a) , 
    
    v = u^2 +ae  ,


for u,v>0. Notice that (<ref>) defines a convex parabola with vertex at (0,ae), and (<ref>) has an asymptote at u = a /b. For values of u< a/b we have that v<0, which is biologically unfeasible. Now, we find that the nullcline (<ref>) can be recasted as a smooth function of u, which is positive  and posses critical points for u>a/b as v^'(u)=0 is held. In so doing, we get that -2a+b(u-Ku^3)=0 must be satisfied. This yields to, as a consequence of the Descartes' rule of signs, existence of two positive roots such that a/b<u_1^⋆<u_2^⋆ as a,b,K>0, where the lower value corresponds to a local minimum and the greater one to a local maximum. Thus, from both expressions in (<ref>), we get that there are at most three positive steady-states. This is schematically illustrated in Fig. <ref> for a system defined later in (<ref>) which is equivalent to (<ref>) in the interior of the first quadrant and is well-defined along the axis v=0. These sketches provide sensitive evidence that two saddle-node bifurcations may occur as two steady-states are created by varying appropriate parameters; see panels (a) and (c), where the nullclines are tangent. 




 §.§ Local stability  of positive steady-states




We now pay attention to the local stability of positive steady-states (u_*,v_*) in (<ref>). The Jacobian matrix of (<ref>) at (u_*,v_*) is


    𝐉 = ([ 2u_*v_*(1+Ku_*^2)^2-b -u_*^2(1+Ku_*^2)v_*^2;                  2u_*                    -1;                       ]) ,

where, by setting 𝒥:=2u_*/[v_*(1+Ku_*^2)^2]., the trace is given by tra(𝐉) = 𝒥-b-1, and the determinant by det(𝐉) = b-(1-.u_*^2/[v_*(1+Ku_*^2)])𝒥.
Notice that the entries  J_12, and  J_22 of (<ref>) are negative, while J_21, is positive. Hence, local stability features depend on whether the first entry J_11= 𝒥-b is greater or less than one as follows from Hartman–Grobman theorem <cit.>. This is summarised as follows: 


Let be 𝒥:=2u_*/[v_*(1+Ku_*^2)^2]..
The local stability of a positive steady-state (u_*,v_*) of (<ref>) is as follows:



  * It is  locally asymptotically stable  if 𝒥< b+1 and b>(1-.u_*^2/[v_*(1+Ku_*^2)])𝒥.

  * It is a repeller  if 𝒥> b+1 and b>(1-.u_*^2/[v_*(1+Ku_*^2)])𝒥.

  * It is a saddle point if b<(1-.u_*^2/[v_*(1+Ku_*^2)])𝒥.



We now analyse the dynamical behaviour in the neighbourhood of the origin.
 
 



 §.§ Dynamics near the origin





 Despite vector field (<ref>) being defined in 𝒟 = {(u,v)∈ℝ^2 |  u ≥ 0,   v >0 }, a solution with v=0 corresponds to absence of static bacteria, which is a feasible  biological scenario. In consequence, similarly to (<ref>) we consider a map t↦  t / v to get

    Y_0:{[             u̇ = f̃(u,v) , f̃(u,v)=u^21+Ku^2+av-buv ,;             v̇ = g̃(u,v) ,     g̃(u,v)=u^2v+aev-v^2 , ].

which is topologically equivalent to (<ref>) and can be continually extended to the boundary axis v=0.
 
 

The origin of (<ref>) is a repeller.



 Since the eigenvalues of the Jacobian matrix at the origin are λ_1=0 and λ_2=ae>0, with associated eigenvectors v_1=(1,0)^t and v_2=(1/e,1)^t, respectively, it follows that there exists a centre manifold


    W^c={(u,v): v=h(u):=a_2u^2+O(u^3), |u|<δ},


for δ>0 sufficiently small <cit.>.  Since W^c is invariant one must have


    g̃(u,h(u))≡ h'(u)f̃(u,h(u)),


for every (u,v)∈ W^c. It follows from this identity that a_2=0 and h(u)=O(u^3). Hence, the dynamics of (<ref>) restricted to W^c is equivalent to ẋ=x^2+O(x^3). Therefore, the origin of (<ref>) is a repeller.






§ BIFURCATION ANALYSIS  IN THE 2D SYSTEM




 We now perform a numerical exploration with MatCont <cit.> in order to identify the most relevant bifurcations when parameters e or b are  varied slowly: Upon varying  b, different decaying bacteria rate scenarios are explored; meanwhile, variations of e reveal the impact of having several response rates as well as different production rates. We focus our analysis on the way bacteria decaying rates coordinate the emergence of oscillations. 

Figure <ref> depicts bifurcation diagrams of system (<ref>) for parameters e and b. In the left-hand panel,  two saddle-node (SN) points are identified at e_1 and e_2 (see Fig. <ref>(a)), where the system posses only one unstable or stable steady-state for e<e_1 and e>e_2, respectively. On the other hand, in Fig. <ref>(b), two SN points are shown at b_1 and b_2 in addition to a Hopf point (denoted by H, for simplicity) at b_3. This Hopf bifurcation  is  supercritical as the first numerically computed Lyapunov coefficient is negative. Hence, for b_1<b<b_3 a stable limit cycle arises  and coexists with the stable equilibrium at the lower branch. This opens the door for the emergence of  synchronising dynamics for b<b_3 under suitable initial conditions.





A two-parameter continuation from the SN points in Fig. <ref>(a) as parameters  e and a are slowly varied reveals a cusp point (CP) and a Bogdanov–Takens point (BT); see Fig. <ref>. 
In addition,  a two-parameter continuation in Fig. <ref> in terms of both b and e starting from the H point shows a Bogdanov–Takens bifurcation as the organising centre of the dynamics. 
Upon taking into account Proposition <ref>, the parameter space of interest for b and e can be divided into the following regions: 

[label=(*)]

  (I) Unstable Origin and Stable Steady-State: there is  a stable steady-state and the origin is unstable; for values of b close to naught the positive steady-state is a stable node and is a stable focus for slightly larger values, which is a consequence of a transition from negative real eigenvalues to complex ones with negative real part. 

  (II) Oscillation Dynamics: The occurrence of positive stable periodic solutions shape synchronising dynamical features in the bacteria population; this feature suggests a coordinated communication among bacteria; 

  (III) Unstable Spiral, Saddle and Stable Node: there are three steady-states from which only one of them is stable. Moreover, notice that region III is subdivided by a bold dashed curve  of homoclinic bifurcation emerging from the BT point:  Limit cycles exist on the left region (adjacent to II); 

  (IV) Bistability: this region is composed of  three steady-states: a stable node, a stable focus and a saddle. 







Notice that Fig. <ref> also shows that the larger the value of e, the smaller the range of values of b that can sustain oscillations. In other words, as e surpasses a certain value no synchronising dynamics are possible, no matter the value of b. Parameter e plays a bacteria production role in the model, it is directly proportional to the fraction of static bacteria production rate and inversely proportional to the constant autoinducer concentration. This could indicate that the oscillatory dynamics of the population is dependent on the interaction between the decaying and production rates of the bacteria or the decaying rates and the autoinducer concentration. The latter scenario will be explored throughout section <ref>.




 §.§ Hysteresis and synchronisation



Now, as we have found such key dynamical events, which are characterised by parameters b and e, we illustrate these transitions by considering a time-dependent varying variable b=b(t). 





  §.§.§ Linearly varying b




The model shows hysteresis as b is varied linearly forth and backwards. For instance, in Fig. <ref>(a) we fix e=1.0 and continuously vary linearly b=1.0 up to b=3.0. As can be seen, the solution decays in oscillatory fashion  for values of b small enough; then, as b crosses a corresponding Hopf bifurcation value, the solution amplitude starts increasing approaching a stable limit cycle; then it finally approaches a stable steady-state after the homoclinic bifurcation. In Fig. <ref>(b), parameter b is traversed backwards, which yields a solution that remains in a stable node steady-state throughout the first two regions it crosses. As it reaches the Oscillation Region, the solution spikes and approaches a limit cycle. Finally, the solution begins to show a decaying amplitude oscillatory behaviour. 

Notice that in Fig. <ref>(b), the solution's behaviour is rather unsimilar to the one when b increases (see Fig. <ref>(a)), which indicates the existence of hysteretic-like features in the system as b is varied forth and backwards. That is, solutions departs and returns to the same initial state visiting different dynamical regimes depending on parameter's pathways.  In other words, as solution show either oscillatory and non-oscillatory properties with distinctive shape features, this observation may be interpreted as a robust quality of such interaction, where this special feature may be  explained biologically as the ability to maintain functionality against external and internal perturbations;  see, for instance, <cit.>. Varying b from b=1.0 to b=3.0 and backwards illustrates how the dynamics of the solution change for different decaying bacteria rate ratios for a constant value of e. When b values are close to one, the decaying rates of constant and static bacteria are similar. As b grows, the difference between the decaying rates of the bacteria increases.

On the other hand, observe that in Fig. <ref>(c), parameter e=1.5  is kept fixed and b varies from b=1.0 up to b=1.5 and remains at this value as the solution converges to a periodic solution. This periodic behaviour shows that the concentration of motile bacteria slightly rises before the static bacteria concentration and also, its exponential decay is faster than the static population one.
 
 




  §.§.§ Periodically varying b




The bifurcation diagram in Fig. <ref> indicates that parameter b plays a key role in the emergence of oscillatory behaviour. We therefore proceed to illustrate crucial consequences as b varies periodically back and forth between b=1.0 and b=3.0 in a given period of time for a fixed e value. In doing so, we will gain insight on the effect of bacteria decay rates that vary in time. We consider an interval given by t  ∈  [0,2000] and a periodic decay rate parameter given by b(t)=A+Bsin(ω t + φ), where A and B determine parameter bounds and ω and φ are the frequency and phase at which b varies. Fig. <ref> shows three scenarios for b varying at different frequencies. These simulations suggest that, in order to observe meaningful dynamics, the value of ω must be small.
For small ω,  the solution will indeed show oscillations with a significantly small amplitude, but as the value of ω increases, the solution will behave as a constant solution. In Fig. <ref> the solutions for ω = 0.010,  0.027 and 0.044 are depicted. For ω= 0.010, we see that as b starts to increase, the solution oscillates regularly. 
 As b(t) varies periodically the solutions alternate between periodic behaviour, oscillatory decay, and constant. That is, as b increases and  the solution enters the oscillation region in Fig. <ref>, the amplitude of the oscillations increases to remain constant. Notice that, as b starts to decrease, the solution becomes constant until it re-enters the oscillation region, when it spikes. This is characterised by the valley-like windows; see upper panel. When ω=0.027 and 0.044 we see that, as parameter b takes decreasing values, the solution spikes as b enters the oscillation region. However, as ω increases, the stable spiral behaviour lasts less time, until it is no longer present in the solution. These sudden spikes in the solution indicate excitatory dynamics in the system. The previous results convey that the dynamics of the 2D-system (<ref>) is sensible to fast varying parameter b, which suggests that the timescale at which the bacteria decay rates change has a rather direct impact on the possible emergence of synchronisation dynamics. 






 §.§ Andronov–Hopf bifurcation in the 2D system




We now proceed to compute the normal forms for the Hopf bifurcation as is crucial for the dynamical shaping for the distinguished parameter values of Table <ref>. 

Let (U,V)∈ℝ^2_+ be the positive coordinates of an equilibrium of (<ref>) in the first quadrant. Then the set of equations f(U,V)=0 and g(U,V)=0
define implicitly a locally invertible transformation given by Ψ:Λ⟶ℝ^4_+,


(U,V,b,K)↦(a,e,b,K):=   (U (b (1 + K U^2) V-U)V + K U^2 V, (V-U^2) (V + K U^2 V)U (b (1 + K U^2) V-U), b,K),


in Λ:={(U,V,b,K)∈ℝ^4_+: b (1 + K U^2) V-U>0,   V-U^2>0}.
Then the vector field (<ref>) in parameter space Λ has the form:


    {
    x'    =	 x^2y + K x^2 y + U (b (1 + K U^2) V-U)V + K U^2 V - b x ,  
    
    y'    =	 x^2 + V -U^2 - y, 
    .


where, for convenience, we use notation (x,y) to name the state variables, and we denote x'=dx/dt, y'=dx/dt. System (<ref>) is C^∞-equivalent to (<ref>) in parameter space  Λ. Moreover, the (positive) equilibrium coordinates appear now as the explicit parameters (U,V). In what follows in this section, we will give conditions such that  (<ref>) undergoes a generic Hopf bifurcation at (U,V); see <cit.> for more details.

The Jacobian matrix of (<ref>)  at (U,V) is


    J(U,V)=[ -b+2 U(1 + K U^2)^2 V   -U^2(1 + K U^2) V^2;                      ;                    2U                    -1 ].


The trace T and determinant D of J(U,V) are given by


    T=T(U,V,b,K)=-1 -b+2 U(1 + K U^2)^2 V


and

D=[J(U,V)]= (V + K U^2 V)^-2D̂, respectively, where

    D̂=D̂(U,V,b,K)=2 U^3 + 2 K U^5 - 2 UV + b V^2 + 2 b K U^2 V^2 + b K^2 U^4 V^2.


Whenever T=0 and D>0, the eigenvalues of J(U,V) are purely imaginary and non-trivial. Moreover, since we have partial derivative T_b(U,V,b,K)=-1<0,

the Hopf bifurcation in (<ref>) is generically unfolded by parameter b. In particular, it follows from there that equation T(U,V,b,K)=0 implicitly defines the function b(U,V,K)=-1 +2 U(1 + K U^2)^2 V.

We now calculate the  first Lyapunov quantity <cit.> in order to determine genericity conditions. We follow
the derivation in <cit.> and move the equilibrium (U,V) of the system (<ref>) to the origin via the translation x↦ x+U, y↦ y+V to obtain the equivalent system


    {
    x'    =	 U (b - UV + K U^2 V) - b (U + x) + (U + x)^2V + y + K (U + x)^2 (V + y),  
    
    y'    =	 2 U x + x^2 - y .
    .


In particular, the Jacobian matrix of (<ref>) at the equilibrium (0,0) coincides with J(U,V) in  (<ref>).
 Upon substituting the function b(U,V,K) into J(U,V), we get
 

    J_H(U,V)=[                   1 -U^2(1 + K U^2) V^2;                    ;                  2U                  -1 ],


with T_H:=[J_H(U,V)]≡0 and [J_H(U,V)]=D_H:=D|_T=0=2 U^3 - V^2 - K U^2 V^2(1 + K U^2) V^2.
 
If D_H>0, then 𝐯_1=(12U,1)^t and 𝐯_2=(-w2U, 0)^t are the generalised eigenvectors of J_H(U,V),
where w=√(D_H).

The change of coordinates (x,y)^t↦ [𝐯_1 𝐯_2] (x, y)^t, where t stands for transpose, allows us to express system (<ref>) with T_H=0 in the form


    ([ x'; y' ])=( [  0 -w;  w  0 ])([ x; y ])+([ P(x,y); Q(x,y) ]),


where



    P(x,y)=14 U^2x^2 - w  2 U^2 x y+ w^2 4 U^2y^2,

and

    [                                                                                                          Q(x,y)                                                                                                               =                   -16 K U^7 - 8 K^2 U^9 - 2 U V^2 + V^3 + 3 K U^2 V^3 + 3 K^2 U^4 V^34 U^2 (V + K U^2 V)^3 wx^2;                                                                                                                                               +  K^3 U^6 V^3 + 8 U^5 (-1 + K V) + 2 U^3 V (4 + 3 K V)4 U^2 (V + K U^2 V)^3 wx^2;                                                                                                                 +-4 K U^5 + 2 U V - V^2 - 3 K U^2 V^2 - 3 K^2 U^4 V^2 - K^3 U^6 V^2 -
     U^3 (4 + 6 K V)2 U^2 (1 + K U^2)^3 V^2xy;                                                                                                                                           +(-2 U + 6 K U^3 + V + 3 K U^2 V + 3 K^2 U^4 V + K^3 U^6 V) w4 U^2 (1 + K U^2)^3 Vy^2;                                                                                                                                        +  12 K^2 U^8 + 4 K^3 U^10 - 4 K U^6 (-3 + K V) + V^2 (1 + 2 K V)2U (V + K U^2 V)^4 wx^3;                                                                                                                                           +    U^4 (4 - 8 K V - 3 K^2 V^2) - 2 U^2 V (2 + K V + K^2 V^2)2U (V + K U^2 V)^4 wx^3;                                                                                                                     + 4 K^2 U^6 - 2 V (1 + 3 K V) + 2 K U^4 (4 + 3 K V) + 
     U^2 (4 + 4 K V + 6 K^2 V^2)2U (1 + K U^2)^4 V^3x^2y;                                                                                                                                                        +(1 - 2 K (U^2 - 3 V) - 3 K^2 (U^4 + 2 U^2 V)) w2U (1 + K U^2)^4 V^2xy^2;                                                                                                                                                                        + K (-1 + K U^2) w^2U (1 + K U^2)^4 Vy^3+O(||(x,y)||^4),;                                                                                                                 ]


 
 in a Taylor expansion near (x,y)=(0,0).
 
 System (<ref>), equations (<ref>) and w=√(D_H) allow us to use the derivation in <cit.> for the direct calculation of the first Lyapunov quantity L_1. In so doing, we obtain the following expression:


    L_1= 164 w^2 (1 + K U^2)^6 U^4V^5  l_1,

where


    [                                                                                                l_1                                                                                                  =                     -32 K^3 U^14 + 16 K^4 U^13 V^2 w^2 + 2 U^2 V^3 (2 + 3 K V^2 (w-1) w) (1 + w^2);                                                                                                                                                  - 2 U V^4 (1 + w^2)^2 +
       V^5 w ( -1+w - w^2 + w^3);                                                                                                       - 
      4 U^3 V^3 (-1 + (1 + 6 K V) w^2 + 6 K V w^4)+ 
      4 K^3 U^11 V^2 (16 w^2 + K (V + 7 V w^2));                                                                                                              + 
      U^4 V^2 (-24 K V (1 + w^2) - 8 (3 + w^2) + 
         15 K^2 V^3 w (-1 + w - w^2 + w^3));                                                                                                    + 
      4 U^6 V (12 + 9 K^2 V^2 (1 + w^2) + 4 K V (3 + w^2) + 
         5 K^3 V^4 w (-1 + w - w^2 + w^3));                                                                                                                                          + 
      K^2 U^12(-96 - 48 K V + K^4 V^5 w (-1 + w - w^2 + w^3));                                                                                                                                           + 
      6 K U^10(-16 - 8 K V + K^4 V^5 w (-1 + w - w^2 + w^3));                                                                                                               + 
      U^8 (-32 + 48 K V + 24 K^2 V^2 (3 + w^2)+ 
         15 K^4 V^5 w (-1 + w - w^2 + w^3));                                                                                                                                      - 
      4 U^5 V^2 (-4 w^2 - 4 K V (1 + w^2) + 3 K^2 V^2 (-1 + w^4));                                                                                                                 + 
      2 K^2 U^9 V^2 (48 w^2 + 8 K (V + 5 V w^2) + 
         3 K^2 V^2 (1 + 6 w^2 + 5 w^4));                                                                                                                     + 
      8 K U^7 V^2 (8 w^2 + 3 K (V + 3 V w^2) + 
         K^2 V^2 (2 + 7 w^2 + 5 w^4)).;                                                                                                    ]


Thus, we have obtained the following result.


Let (U,V,b,K)∈Λ be such that T_H=0, D_H>0 and l_1≠0.
Then (<ref>) undergoes a codimension-one Hopf bifurcation at the equilibrium (U,V). In particular, if l_1<0 (resp. l_1>0), the Hopf bifurcation is supercritical (resp. subcritical), and a stable (resp. unstable) limit cycle bifurcates from (U,V) under suitable parameter variation. 


Whenever condition l_1≠0 in Theorem <ref> does not hold, the Hopf bifurcation of (<ref>) at (U,V) is degenerate. The actual codimension of this singularity —and the stability of further limit cycles that bifurcate— is determined by the sign of the so-called second Lyapunov quantity L_2; see, for instance, <cit.>. However, since the transformation (<ref>) is invertible in the parameter set Λ, numerical evidence suggests that l_1<0 for representative parameter values in Table <ref> and, hence, the bifurcating limit cycle is stable.
 Furthermore, bifurcation theory ensures that the existence and stability of this stable periodic orbit persist for an open set of parameter values in the region T(U,V,b,K)>0, i.e., when the focus at (U,V) is unstable <cit.>.





 §.§ Bogdanov–Takens bifurcation in the 2D system



We give now conditions such that our model undergoes a Bogdanov–Takens bifurcation under suitable parameter variation at a positive equilibrium. We prove the existence of a germ of a BT bifurcation and show that our system (under certain conditions) is locally topologically equi­va­lent to a normal form of the BT bifurcation. We refer to <cit.> and the references therein for the derivation of the genericity and transversality conditions that need to be verified during this proof.

For the sake of clarity, it is convenient to state the dependence of the vector field (<ref>) on parameters b and K explicitly. Hence, throughout this section we denote X: ℝ^4_+⟶ℝ^2_+,


    X(x,y;b,K) =  ( x^2(1 + K x^2)y + α - b x, 
     x^2 + αϵ - y  ),


where we use notation (x,y) for the state variables.
Also, let us denote the Jacobian matrix of X with respect to the variables (x,y) as

    ∂ X∂(x,y)(x,y;b,K)=(
    [ --2 x + b y + 2 b K x^2 y + b K^2 x^4 y(1 + K x^2)^2 y                                    -x^2(1 + K x^2) y^2;                                                     2x                                                     -1;                                                        ]
    )
    .



Step 1. We verify that the system may exhibit a singularity with a double zero eigenvalue and geometric multiplicity one. 


Consider the set

    [                                                                    Ω_BT                                                                       =   {(x,y)∈ℝ^2_+:  -2 x^5 + y^3>0, 2 x^3 - y^2>0, -2 x^5 + y^3- x^3 y >0,;                                                                         (x^2 - y)(2 x^5 + x^3 y - y^3)>0, 4 x^5 - x^3 y - 6 x^2 y^2 + 4 y^3≠0}. ]


Then there is a positive equilibrium (U,V)∈Ω_BT of  (<ref>) and there are positive parameter values (b^*,K^*)=(-2 U^5 + V^32 U^5,2 U^3 - V^2U^2 V^2)∈ℝ^2_+ such that .∂ X∂(x,y)|_(U,V;b^*,K^*) is nilpotent. 




In order to prove this Lemma, we look look for solutions (x,y,b,K) of the algebraic system



    x^2(1 + K x^2)y + α - b x   =   0 ,
    
    x^2 + αε - y    =   0,
    [∂ X∂(x,y)(x,y;b,K)]   =   0,
    [∂ X∂(x,y)(x,y;b,K)]   =   0,



with (x,y)∈Ω_BT, b>0 and K>0.
Solving for b and K in (<ref>)-(<ref>) leads to 
 b=-2 x^5 + y^32 x^5>0 and K=2 x^3 - y^2x^2 y^2>0. 
 Substitution of  (b,K) in (<ref>)-(<ref>) and solving for (α,ε) leads to
 

    α=-2 x^5 + y^3- x^3 y 2 x^4  and  ε =2 x^4 (x^2 - y)2 x^5 + x^3 y - y^3 .


In particular, notice that  α>0 and ε>0 since (x,y)∈Ω_BT.  The map (x,y)↦(α,ε) defined by (<ref>) has a Jacobian matrix given by 


    ∂(α,ε)∂(x,y)=(
    [                                                          -1 + y2 x^2 - 2 y^3x^5                                                               -x^3 - 3 y^22 x^4; 2 x^3 (2 x^7 + 5 x^5 y - x^3 y^2 - 6 x^2 y^3 + 4 y^4)(2 x^5 + 
      x^3 y - y^3)^2                           -6 x^9 + 6 x^6 y^2 - 4 x^4 y^3(2 x^5 + x^3 y - y^3)^2;                                                                                 ]
    ) ,

which determinant is

    [∂(α,ε)∂(x,y)]=4 x^5 - x^3 y - 6 x^2 y^2 + 4 y^3x(2 x^5 + x^3 y -  y^3).

Therefore, since [ ∂(α,ε)∂(x,y)]≠0, the Inverse Function theorem ensures that system (<ref>) is locally invertible. Then, the solution of (<ref>) is given by (x,y,b,K)=(U,V,b^*,K^*)∈ℝ^4_+, where (U,V) is locally defined by the inverse map of (<ref>) and b^*=-2 U^5 + V^32 U^5 and K^*=2 U^3 - V^2U^2 V^2.

 
Hence, at (b,K)=(b^*,K^*) the equilibrium (U,V) of the system (<ref>) has a Jacobian matrix given by


    ∂X∂(x,y)(U,V;b^*,K^*)=(
    [    1 -12U;   2U   -1;      ]
    )

with a double zero eigenvalue. In particular, note that (<ref>) is not the null matrix. 
The corresponding generalized eigenvectors of (<ref>) are given by 


    𝐯_1=(12U,1)^t, and𝐯_2=(1,-1 + 2 U)^t .


It follows that (<ref>) is nilpotent and that the double zero eigenvalue has geometric multiplicity one.

 


Step 2.
The next goal is to state the following transversality condition of a Bogdanov–Takens bifurcation.

 Let (U,V,b^*,K^*) be as in Lemma <ref> and 
consider the map Ψ:ℝ^4→ℝ^4,

    (x,y,b,K)↦(x^2(1 + K x^2)y + α - b x, 
     x^2 + αϵ - y  ,T,D),

where T and D are the trace and determinant of the matrix 

∂ X∂(x,y)(x,y;b,K),

respectively. Then the map Ψ is regular at (x,y,b,K)=(U,V,b^∗,K^∗).



The 4×4 Jacobian matrix DΨ=DΨ(x,y;b,K) of the map Ψ is


    DΨ=(
    [               -b - 2 K x^3(1 + K x^2)^2 y + 2 x(1 + K x^2) y                                          -x^2(1 + K x^2) y^2                                                           -x                                        -x^4(1 + K x^2)^2 y^2;                                                           2x                                                           -1                                                            0                                                            0;                                   2 - 6 K x^2(1 + K x^2)^3 y                                          -2 x(y + K x^2 y)^2                                                           -1                                        -4 x^3(1 + K x^2)^3 y; 8 K x^4 + 2 K^2 x^6 - 2 y + 6 x^2 (1 + K y)(1 + K x^2)^3 y^2                    -4 x^3 - 4 K x^5 + 2 x y(1 + K x^2)^2 y^3                                                            1                  -2 (x^5 + K x^7 - 2 x^3 y)(1 + K x^2)^3 y^2;                                                              ]
    ).


After some calculations, we have [DΨ(x,y,b,K)]=2x^5(y + K x^2 y)^4F(x,y,b,K) with

    F(x,y,b,K)=-2 K x^5 - 9 x y + b y^2 + 2 b K x^2 y^2 + b K^2 x^4 y^2 + 
     x^3 (10 + K y).
 
In particular, straightforward substitution and algebraic simplification leads to 
F(U,V,b^∗,K^∗)=-2UV^2(4 U^5 - U^3 V - 6 U^2 V^2 + 4 V^3)≠0, since (U,V)∈Ω_BT.
It follows that [DΨ(U,V,b^∗,K^∗)]≠0,
which ensures that the map Ψ is regular at (x,y,C,Q)=(U,V,b^∗,K^∗).

 

Step 3. 
We now construct a change of coordinates to transform X(x,y;b,K) into a normal form of the Bogdanov–Takens bifurcation; we refer to <cit.> once again.

 Let (U,V,b^*,K^*) be as in Lemma <ref> and consider the following auxiliary expressions:


    G_1=8 U^10 - 2 U^8 V - 4 U^5 V^3 - 3 U^3 V^4 + 2 V^6 ,    G_2=2 U^5 + 3 U^3 V - 2 V^3 .


If G_1≠0 and G_2≠0, then there exists a smooth, invertible transformation of coordinates, an orien­tation-preserving time rescaling, and a re­pa­ra­me­tri­za­tion such that, in a sufficiently small neighbourhood of (x,y,b,K)=(U,V,b^∗,K^∗),  system (<ref>) is topologically equivalent to a normal form of the codimension-two Bogdanov–Takens bifurcation.



Let us set up the equilibrium (U,V)∈Ω_BT of (<ref>) to the origin via the translation x↦ x+U, y↦ y+V to obtain the equivalent system


    Y:{
    x'    =	 α  + (x + U)( -b+ (x + U)(1 + K (x + U)^2) (y + U)),  
    
    y'    =	 αε + (x + U)^2 - y - V .
    .


In particular, the Jacobian matrix of (<ref>) at the equilibrium (0,0) at the bifurcation point (b^*,K^*) coincides with ∂X∂(x,y)(U,V;b^*,K^*) in  (<ref>).

  Let 𝐏=[𝐯_1, 𝐯_2] be the matrix whose columns are 𝐯_1 and 𝐯_2; see (<ref>). Next, consider the following change of coordinates:


    (
    [ u; v;   ]
    )
    =
    𝐏^-1(
    [ x; y;   ]
    ).


Then, the vector field given by


    𝐉=𝐏^-1∘ Y∘𝐏,


is 𝒞^∞-conjugated to Y in (<ref>). 

Taking a Taylor expansion of 𝐉(u,v;b,K) with res­pect to (u,v) around (u,v)=(0,0) and evaluating at (b,K)=(b^∗,K^∗), one obtains


    (
    [ u̇; v̇;    ]
    )
    =
    (
    [ 0 1; 0 0;   ]
    )
    (
    [ u; v;   ]
    )+
    12(
    [ a_20u^2+2 a_11uv+a_02v^2+O(||(u,v)||^3); b_20u^2+2 b_11uv+b_02v^2+O(||(u,v)||^3);                                         ]
    ),


where we have that: a_20=14U^10 V( 8 U^10 - 16 U^11 + 4 U^9 V - 4 U^5 V^3 + 8 U^6 V^3- 3 U^3 V^4 + 6 U^4 V^4 + 2 V^6 - 4 U V^6), b_20=14 U^10 V(8 U^10 - 2 U^8 V - 4 U^5 V^3 - 3 U^3 V^4 + 2 V^6) and b_11= 12 U^9 V(-4 U^9 + 8 U^10 - 2 U^8 V + U^4 V^3) - 12 U^9 V(4 U^5 V^3 + 3 U^3 V^4 - 2 V^6).


If b_20≠0 and a_20+b_11≠0, then the theory of normal forms for bifurcations <cit.> ensures that our system fulfills the necessary genericity conditions to undergo a codimension two Bogdanov–Takens bifurcation. 
In particular, condition G_1≠0 ensures that b_20≠0. Furthermore, after some algebraic manipulation one obtains a_20+b_11= -V^2G_2/4U^10.

In summary, Lemmas <ref> and <ref>, and inequality G_1G_2≠0  ensure that the genericity and transversality conditions of a codimension two Bogdanov–Takens normal form are satisfied. Hence there exists a smooth, invertible transformation of coordinates, an orientation-preserving time rescaling, and a reparametrization such that, in a sufficiently small neighbourhood of (x,y,α,k)=(U,V,b^∗,K^∗), the system (<ref>) is topologically equivalent to one of the following normal forms of a Bogdanov–Takens bifurcation:


    {
    [                       ξ̇_1                          =                       ξ_2,;                       ξ̇_2                          = β_1+β_2 ξ_2 +ξ_2^2±ξ_1ξ_2,;                            ]
    .


where the sign of the term ξ_1ξ_2 in (<ref>) is determined by the sign of (a_20+b_11)b_20.




 



The next theorem is a straightforward consequence from the findings in Lemmas <ref>, <ref> and <ref>, as can be seen in <cit.>, and it summarises the main result in this section.



Let (U,V,b^*,K^*) be as in Lemma <ref> and consider the quantities G_1 and G_2 defined in Lemma <ref>.
	Then if (b,K)=(b^*,K^*), system (<ref>) undergoes a codimension-two Bogdanov–Takens bifurcation at (x,y)=(U,V).







§ BIFURCATION AND CHAOS IN THE 3D SYSTEM




In order to understand the global picture that emerges, we perform a bifurcation analysis of system (<ref>) with Auto <cit.> allowing both parameters k_2 and γ to slowly vary. The other parameters remain fixed in their typical values, except for μ_1=0.2, and μ_2=0.7. Fig. <ref> shows the bifurcation diagram in the (k_2,γ)-plane. There curves of saddle-node (LP) and supercritical Hopf (labelled in this section as HB) bifurcation meet at a BT point. Also, a curve of Shilnikov homoclinic bifurcation (hom) emerges from the point BT. A codimension-two Belyakov point (B_c) on the curve hom marks the onset of chaotic dynamics. To the right of B_c, the homoclinic bifurcation is simple; and to the left of B_c, the homoclinic bifurcation is chaotic. More concretely, one can find horseshoe dynamics in return maps defined in a neighbourhood of the homoclinic orbit. The suspension of the Smale horsehoes form a hyperbolic invariant chaotic set which contains countably many periodic orbits of saddle-type.
The horseshoe dynamics is robust under small parameter perturbations; hence, the chaotic dynamics persist if the homoclinic connection is broken; see <cit.>.
A second codimension-two point B_+ and also called a Belyakov point, lies on hom for smaller values of both k_2 and γ. At the point B_+, the steady-state associated with the homoclinic orbit has repeated stable eigenvalues.
On the  segment of hom to the right of  B_+, the homoclinic orbit converges to a saddle-focus (i.e., it has a complex pair of stable eigenvalues); and on the segment to the left of B_+, the same equilibrium is a real saddle (i.e., the stable eigenvalues are real). 

 



The bifurcation picture in Fig. <ref> is just a partial representation of the full complexity one may encounter in this region of parameter space. Indeed, the saddle periodic orbits may also undergo further bifurcations such as period-doubling and torus bifurcations <cit.>. Moreover, the presence of the chaotic Shilnikov homoclinic bifurcation and that of the Belyakov points B_c and B_+ imply a very complicated structure (not shown) of infinitely many saddle-node and period-doubling bifurcations of periodic orbits as well as of subsidiary n-homoclinic orbits <cit.>. Moreover, for each of these subsidiary n-homoclinic orbits, the system exhibits countably many horseshoes as in the original homoclinic scenario.  


The Belyakov points B_c and B_+ in Fig. <ref> divide the curve hom in three segments. We now fix parameters (k_2,γ) at selected points on each of these segments —labelled as ♢, △, and ∘, respectively— and allow parameters μ_1 and μ_2 to vary. The resulting picture is shown in Fig. <ref> in the top row, while suitable enlargements are presented in the bottom row.  While  the resulting bifurcation scenarios shown in Fig. <ref> are slightly different to one another, the main ingredients organizing complicated dynamics (such as Belyakov points and homoclinic chaos) remain a common feature in each case. In particular, in the three cases, the bifurcation curves corresponding to the scenario of Fig. <ref> are those occurring for lower values of μ_2 in Fig. <ref>. The second curve LP  (that which is present for larger values of μ_2) and associated bifurcation phenomena is not present in Fig. <ref>.
The Hopf bifurcation curve hom in Fig. <ref> is now separated into two segments by a degenerate Hopf point GH from which a curve of saddle-node of periodic orbits (LPC) emerges.
The curve of homoclinic bifurcation hom now contains further codimension two resonant (R) and Belyakov (B_-) points. While a single stable cycle bifurcates from B_- (and no extra bifurcations occur),  the full picture near the points R include period-doubling bifurcations and a curve of 2-homoclinic bifurcation (not shown) emanating from the codimension-two points. Of particular interest is the case in panels (c) where the curve hom terminates at a BT point located on the second curve LP.  Also, from this BT point a  subcritical Hopf bifurcation emerges. In particular, this HB curve is very close to the saddle-node curve; these two curves meet at a Zero-Hopf bifurcation point (ZH). 
The exact dynamical features which appear in phase space for parameter values near a ZH point depend on higher order terms in a normal form approach, which is beyond the scope of this work. It suffices to say that this codimension-two point is often related to the emergence of invariant tori and chaotic invariant sets <cit.>. 
For larger values of μ_2 and lower values of μ_1, another (non chaotic) homoclinic bifurcation curve makes a sharp turn near the point ZH.


 




Finally, Fig. <ref> shows bifurcation diagrams in the (k_2,γ)-plane with values of μ_2=1 (in panel (a)), and μ_2=1.4 (in panel (b)). The resulting bifurcation diagrams are similar to that in Fig. <ref>, but with higher values of μ_2. In  Fig. <ref> the HB curve has a degenerate Hopf point GH, from which a curve of saddle-node of cycles LPC emerges. All in all, upon comparing figures <ref> and <ref>, an increase of parameter μ_2 produces extra bifurcations that favors chaotic behaviour. Indeed, the hom curve in  Fig. <ref> is now chaotic along its entire length. Also, in Fig. <ref>(b), an increase of μ_2 favors the emergence of a ZH point. However, the bifurcation sets are “pushed” towards lower values of k_2 and larger values of γ as μ_2 is increased; compare the scale of the variables in figures <ref> and <ref>. 





§ DISCUSSION



Analysing the effect of autoinducers on bacterial populations is crucial due to their key role in the QS mechanism. We presented a model that demonstrates the interaction between autoinducers and two subtypes of bacteria. A thorough analysis of this model revealed parameter sets that lead to oscillation dynamics in both the constant autoinducer sub-model and the full three-component model.
The study was performed with a combination of both rigorous normal form analysis and numerical continuation methods.
 Upon analysing the full system, we were able to understand how autoinducer concentrations interact with bacterial populations and influence bacterial communication triggered by these interactions.

The constant autoinducer system (<ref>) shows that parameter e is inversely proportional to the autoinducer concentration q_0, indicating that changes in e can impact response rates and production rates on the population. Exploring additional parameters revealed that parameter b, related to bacteria decaying rates, leads to a family of limit cycles via Hopf bifurcation. By conducting two-parameter numerical continuation on parameters b and e, the (b,e)-plane was divided into four regions based on their dynamical events. It was observed that no synchronising behaviour occurs for high values of e, suggesting that oscillatory behaviour may not occur below certain critical shallow threshold autoinducer concentration; this discovery agrees with QS features already recognised in e. g. <cit.>.

Furthermore, as can be seen from Proposition <ref>, when 𝒥=b/(1-. u_*^2/[v_*(1+Ku_*^2)])., notice that tra(𝐉)=[(b+1)u_*^2 - v_*(1+Ku_*^2)]/[v_*(1+Ku_*^2)-u_*^2]. and det(J)=0, which suggest that other nonlinear events may trigger synchronising behaviour  on the bacterial population not only by means of a Hopf bifurcation. To fully understand such an implication and, in consequence, the impact of autoinducer dynamics on the population, we analysed the 3D system and identified k_2 and γ as key parameters for synchronising dynamics. Through Hopf bifurcations, self-sustained oscillations were observed to emerge. A two-parameter continuation analysis revealed that these parameters play a role in the emergence of additional oscillatory behaviour and robust synchronisation properties in the population.

From section <ref>, it was observed that a Shilnikov homoclinic bifurcation branch originates from a Bogdanov–Takens point, leading to a Belyakov point B_c that separates simple and chaotic behaviour along the homoclinic curve. Another Belyakov point B_+ marks the boundary between real saddle and saddle-focus steady-states. By conducting a two-parameter continuation along three pivotal points on the homoclinic curve, varying parameters μ_1 and μ_2 (equivalent to parameter e in the constant autoinducer system), intricate dynamics were revealed resulting in resonant and additional Belyakov points. Additionally, a zero-Hopf point appears in an Hopf curve intersecting a nearby saddle-node curve. These findings support chaotic dynamics as infinitely many saddle-node, period-doubling orbits and n-homoclinic orbits —accompanied with countably many horseshoes for all n∈ℕ— come to birth as well as invariant tori and invariant chaotic sets.   

The QS model examined in this study depicts synchronised oscillatory dynamics in bacterial populations, potentially exhibiting synchronising characteristics. However, to fully understand the global picture of the QS mechanism, migration dynamics and time delays should also be considered. The inclusion of a time delay is crucial as it influences the dynamics of autoinducer reception and emission, indicating that bacteria do not react instantaneously to autoinducers (e.g. <cit.>). These aspects will be addressed in future research.




§ ACKNOWLEDGEMENTS
  PA thanks Proyecto Interno UTFSM PI-LIR-24-04 and Proyecto Basal CMM-Universidad de Chile. VFBM thanks the financial support by Asociación Mexicana de Cultura A.C.



siamplain



