1 Introduction

Bacterial populations provide an interesting subject for the study of emergence and complexity phenomena. Bacteria react to the environment in complex ways and use short-distance local communication to start coordinated behaviour that can involve hundreds of thousands of cells. Tissues in multicellular organisms, on the other hand, are often controlled centrally via hormones and other signals. The Gram-negative soil-dwelling bacterium, Myxobacteria, best known for population level rippling and fruiting body stages whilst under starving conditions, provides an interesting and extensively studied vehicle for understanding cellular communication and spatial behaviour in bacteria [1,2,3,4,5,6,7]. In particular, phenomena such as swarming, aggregation and rippling in the Myxobacteria Myxococcus xanthus continue to elude researchers in their understanding of the mechanisms by which Myxobacteria communicate and organise themselves to predate prey bacteria [8].

Myxobacteria are able to move over a substrate using slow gliding motions driven by two modes of regulated transport, adventurous (A-motility) and social (S-motility) [7, 9, 10]. Bacteria at the colony edge or under isolation can undergo an adventurous mode of motility through the secretion of slime to propel the organism forward. Through cell-to-cell contact-signal (C-signal) transduction, bacterium undergoes reversal of the intracellular motility mechanism [11]. After a number of hours, the accumulation of C-signals across a population will cause an oscillatory motion within the population, giving a visual impression of a diffusing wave. Although it appears as though colliding waves pass through one another, the head-to-head impact of bacterium causes the wave to reverse resulting in very little net displacement. Under S-motility the cell relies on the extension and retraction of type IV pili: these extend from the leading cell pole and grip the polysaccharide cell surface or solid substrate surface. Both type IV pili-based motility and slime extrusion- based motility require the Frz pathway to induce forward movement in the direction of the long axis of the cell.

P systems are distributed parallel computability models based on cellular membrane compartmentalisation and diffusing chemical signals within and across membranes [12]. A system is defined within a unique outer membrane known as the skin that can hold any number of membrane structures, similar to the organelles inside a living cell. Each membrane has a cellular region which can hold a multiset of objects defined by some alphabet and a set of rewrite rules that enables the multiset of objects to evolve under the correct conditions. Objects can diffuse across membranes analogous to cellular signals, and membranes can dissolve passing their chemical contents onto the parent membrane region. Rules can be assigned a priority rating, giving the set of rules for a single membrane a specific order of execution similar to the downstream transduction of receptor activation. If a high-priority rule cannot execute under current conditions, then any additional rules are tested. The priority of rules is typically denoted as \(r_{1} > r_{2}\), to indicate that the first rule has a higher execution priority over the second rule. Additionally, the prioritising of rules can be assigned using a stochastic constant [13]. These stochastic P systems use the Gillespie algorithm to determine which rule to execute and the order.

The evolution of a P system is performed in parallel across all cellular regions similar to the parallel operation of the organelles within the cell and the chemical operations within the main cytoplasmic space itself. Operations are performed either according to the order of each rule, a stochastic selection process, or a mixture of the two. P systems can generally be divided into two categories: the first, non-cooperative, where only single objects can evolve; and the second, cooperative, where multiple objects can evolve in parallel [14]. If every membrane is incapable of evolving over its multiset, the system is said to halt. Upon halting, values can be read from a specific membrane region as defined in the P system formal definition. Under such conditions, the P system is said to be computationally complete, yet unlike the model proposed in this study, these systems would be used to generate a value or a decision rather than spatial behaviour.

Due to inherent compartmentalisation of a P systems model and their encapsulation of chemical diffusion across membranes, they have proved very apt in modelling Quorum sensing [15,16,17,18,19]. The spatial dynamics of diffusing signals can be captured in a grid-like P system, known as a multi-environment, similar to that of a cellular automaton [20]. Further works of modelling biological systems using membrane computing can be seen in Cardona et al. [21, 22], which look at the dynamics of ecosystems. To address the inherent randomness and uncertainty in large biological systems, multi-environment P systems have been combined with probabilistic terms [23, 24].

We explore the use of the P system formalism to increase our understanding of spatial phenomena in Myxobacteria. Recent extensions to the P system formalism allow simulating the models with a stochastic engine, which is crucial to capture realistic biochemical scenarios [25]. While the spatial behaviour of Myxobacteria is controlled primarily by a non-diffusible contact-based signal, we explore commonalities between this form of communication and chemical quorum sensing which underlies other modelling work using P systems.

2 Methods

The methods proposed in this study begin with a brief description of the formal definition for a P system; the reader is advised to refer to the work of Păun et al. [26]. The P system construct is then described in light of an extension to incorporate a two-dimensional spatial construct. The execution of such a multi-environment P system is covered with a brief description of the Gillespie algorithm. Finally, the two-dimensional model and Myxobacterial rewrite rules are reviewed.

2.1 P system formal definition

P systems can be graphically represented as Venn diagrams and trees. Venn diagrams help illustrate the membrane structure in respects to the physiology of a real biological cell, in addition to showing the multiset of objects and membrane rules.

A P system containing \(n\) cellular regions delimited by membranes is defined formally as

$${{ \Phi}}_{n} = \left( {{{\varSigma }},\mu ,w_{1} , \ldots ,w_{n} ,\left( {R_{1} ,p_{1} } \right), \ldots ,\left( {R_{n} ,p_{n} } \right),i_{0} } \right),$$

where

  • \({{\varSigma }}\) is an alphabet of objects that populate a multiset;

  • \(\mu\) defines the membrane structures including the skin. The labels for the membranes are one-to-one with the number of defined regions. For example, \(\left[ {\left[ {} \right]_{2} \left[ {\left[ {} \right]_{4} } \right]_{3} } \right]_{1}\) defines a P system \({{ \Phi}}_{4}\);

  • \(w_{i}\) defines the initial multiset content of each membrane, as indicated by the respective label. The multiset for each region is defined as the set \(w_{i} \in {{\varSigma }}^{*}\);

  • \(\left( {R_{i} ,p_{i} } \right)\) refers to the rules and rule priority, where \(0 \le p_{i } \le 1\), for the membrane region i. If there is only one rule for a given membrane region the priority is set to 0. Multiple rules for the same region require that each rule be given a priority;

  • \(i_{0}\) defines the output region which should be read only when the system halts.

2.2 Multi-environment P system

A multi-environment P system represents a spatial structure using a grid of membrane regions connected by links. It has the same eight directional properties as a regular two-dimensional cellular automaton using a Moore neighbourhood. The formal definition as presented by Romero-Campero et al. [19]:

$${\text{multi}}\_\text{{environment}} = \left( {{{\varSigma }},H,G,E_{1,1} , \ldots ,E_{x,y} ,R_{1,1} , \ldots ,R_{x,y} ,{{ \Phi}}} \right),$$

which we present as a Myxobacterial environment:

  • \({{\varSigma }} = \left( {{\text{slime}}, \ldots } \right)\) is an alphabet used to construct the contents of a multi-environment region. The alphabet contains chemical signals such as slime and lysing antibiotic chemicals;

  • \({\text{H}}\) is a finite set of labels used to distinguish each multi-environment region. We use the Cartesian ordered pair \(X \times Y = \left\{ {\left( {i,j} \right):i \in X \,\,{\text{and }}\,\,j \in Y} \right\}\) to layout each multi-environment region similar to the coordinates of a two-dimensional cellular automata;

  • \({\text{G}} = \left( {V,S} \right)\) is a graph with nodes \(V = \left( {\left( {1,1} \right), \ldots ,(x,y)} \right)\) represent the multi-environment regions and the edges S define the links between regions. Mycobacteria P systems are able to traverse the multi-environment P system by the link. This is an alternative notation to the P system \(\mu\), which is used to outline a membrane structure;

  • \(E_{i,j} = \left( {w_{i,j} ,R_{i,j} ,\varphi_{n} } \right)\), for each environment region \(\left( {1,1} \right) \le \left( {i,j} \right) \le \left( {x,y} \right)\), \(E_{i,j}\) is the initial configuration of the region. \(w_{i,j} \in {{\varSigma }}^{*}\) is a multiset of characters over the alphabet \({{\varSigma }}\) and \(R_{i,j}\) defines a set of rules and \(\left\{ {\varphi_{n} \in {{ \Phi}}:1 \le n \le 8} \right\}\) defines a number of bacterial P systems within that multi-environment P system cell, where each bacterial P system can adopt up to one of the eight possible directions according to the underlying Moore neighbourhood arrangement of the multi-environment grid structure;

  • \(R_{i,j}\) is a finite set of rule for each multi-environment regions. Rule priority is defined on a stochastic basis, using the likelihood that a rule will be executed given certain environmental constraints;

  • \({{ \Phi}}\) is the inclusion of n Myxobacterial P systems within this particular multi-environment region. A single Myxobacterial P system is denoted by one of eight directions defined by a Moore neighbourhood. The sum of Myxobacteria present in a multi-environment P system is defined by the sum of k (see below) across the eight possible Myxobacterial P systems present at any one time.

2.3 Simulating P system models

The chemical master equation gives the probability of obtaining a state S at time t from an initial set of conditions [25]. Yet, even in short simulations, this form of stochastic modelling can be very computationally expensive [27]. An alternative means of modelling stochastic processes introduced by Gillespie [28] updates the state of a system using probability to determine the order and time for rule execution. The collections of trajectories are updated after the execution of each rule to reflect the change to any of the system variables. In this work, the Gillespie algorithm is modified to calculate the likelihood of each rule being executed in parallel on a constant time step (a unitless value \(\Delta t = 1\)). If a region of the multi-environment P system is not populated with a Myxobacterial P system, the stochastic values for that region are set to zero. The algorithm provides a simple yet elegant means of using environmental properties, such as the volume of a biological space and the probability of chemical collision.

2.4 Myxobacterial P system

The structure of a Myxobacterial P system can be defined as

$$\varphi = \left( {{{\varGamma }},\mu ,L,E} \right),$$

where

  • \({{\varGamma }}\) is an alphabet of multiset elements within the P system region;

  • \({{\mu }}\) represents the structure of a bacterium P system, that is, \(\left[ {\quad} \right]_{l}\) where \(l \in L\);

  • \(L\) is the set \(\left\{ {\text{myxo}} \right\}\) of labels used on the membrane;

  • \(E\) represents the initial configuration of the Myxobacterial P system. Each Myxobacterial P system contains a direction, a means of motility and a bacteria population, defined as n, m and k for the purpose of rule presented in this work.

The configuration of a single Myxobacterial P system is defined by a tuple:

$${{\varGamma ::}} < N,D,M > ,$$

where

  • \(N\), where \(1 \le N \le M\), defines the quantity of Myxobacteria in a P system and M denoted the maximum number of Myxobacterial cells;

  • \(D = \left\{ {{\text{N,}}\;{\text{NE,}}\;{\text{E,}}\;{\text{SE,}}\;{\text{S,}}\;{\text{SW,}}\;{\text{W,}}\;{\text{NW}}} \right\}\) specifies the direction in which a P system of bacteria faces;

  • \(M = \left\{ {A_{\text{m}} ,\;S_{\text{m}} } \right\}\) specifies the motility engine of the cell.

As with the structure of cellular compartmentalisation of a regular P system, the use of multi-environment model in this study provide cellular space for the inclusion of a Myxobacterial P system and the connection between multi-environment regions allows the diffusion of Myxobacterial P systems between multi-environment regions. For example,

$$\left[ {[17,E, - ]_{\text{myxo}} } \right]_{i,j} -\left[ {\quad} \right]_{i + 1,j} \mathop \to \limits^{{c_{i} }}\left[ {\quad} \right]_{i,j}^{\prime } - \left[ {\left[ {17,E, - } \right]_{\text{myxo}} } \right]_{i + 1,j}^{\prime } ,$$

signifies the movement (either A-motility or S-motility, as indicated by the wildcard \(-\)) of a Myxobacteria P system (with a population of seventeen bacterium and orientated toward \(E\)) from a multi-environment region at location \(i,j\), into the adjacent multi-environment region \(i + 1,j\). The superscript \(c_{i} ,\) where \(0 \le c_{i } \le 1\), on the rule transformation refers to the stochastic likelihood of this rule executing. A primed region indicates a change in cellular configuration having executed a rewrite rule.

3 Results and discussion

We begin by describing the behaviour of Myxobacteria as a set of multi-component P system rewrite rules. These rules are then explored further via a demonstration of typical motility, complete S-motility, and complete A-motility from the initial conditions of a bacteria colony. Finally, we see in closer detail the behaviour of Myxobacteria in the presence of slime upon the substrate.

3.1 Stochastic rules

The set of P system rewrite rules are projected over the cellular region of the multi-environment P system and not the Myxobacterial P systems themselves (see Table 1 for the complete list). The list of multi-environment P system rules given here acts as examples of possible cellular states that rewrite rules can operate over the included Myxobacterial P systems. Cellular components of Myxobacteria can change, such as direction and motility, leading to additional emergent behaviour. Rule priority is first determined by whether a Myxobacterial P system is in a suitable state or position. Then a random value [0,1] is calculated and compared to the position of available rules when respective rule priorities are combined and normalised between [0,1]. The availability of some arbitrary n rules implies that \({{\varSigma }}_{i = 1}^{n} c_{i} = 1\).

Table 1 The list of multi-component P system rules
Full size table

The movement of the Myxobacteria is presented in rules \(r_{1}\) through to \(r_{6}\): \(r_{1}\) uses social motility to move a population of Myxobacteria into an empty neighbouring cell; \(r_{2}\) applies to Myxobacteria using adventurous motility to propel itself forward into an empty neighbour whilst depositing slime; \(r_{{ 3 {\text{a}}}}\), \(r_{{ 3 {\text{b}}}}\) and \(r_{{ 3 {\text{c}}}}\) allows the population of Myxobacteria to detect slime perpendicular to its long axis, followed by movement into the slime whilst itself depositing slime; \(r_{{ 4 {\text{a}}}}\), \(r_{{ 4 {\text{b}}}}\) and \(r_{{ 4 {\text{c}}}}\) allows slime detection on an angle; and \(r_{5}\) allows the mixing of a population of Myxobacteria from one neighbouring cell into another, in either motility modes.

The movement of Myxobacteria from one region into another populated region is determined by social motility movement and cell number mixing, perpendicular or on an angle to the moving P system in question, is controlled using \(r_{{ 6 {\text{a}}}}\), \(r_{{ 6 {\text{b}}}}\) and \(r_{{ 6 {\text{c}}}}\).

Collision between populations of Myxobacteria can result in a switch in motility mechanism, from adventurous into social as per C-signal induced reversal (\(r_{7}\)).

As the quantity of adventurous Myxobacteria occupy a single region increases there is an ever-greater likelihood of becoming socially mobile, and therefore joining any social Myxobacteria within that region (as per \(r_{8}\)). The probability of switching motility takes on a rate similar in form to a ligand associating with a substrate as per the Michaelis–Menten equation, and in which a surplus of Myxobacterial asymptotes the probability of switching. On the other hand, if the region is sparely populated with social Myxobacteria, they can switch back into adventurous motility, as seen in \(r_{9}\). Myxobacterium can undergo spontaneous reversal via a switch to the internal cellular motility mechanism (\(r_{10}\)). This causes the Myxobacteria to spend a lot of its solitary, swarming and aggregating life moving along its long axis.

A simplified multi-component P system at time \(t\), as presented in Fig. 1, illustrates some of the rules used to evolve into the configuration at \(t + 1\). Each square represents a multi-environment region whose edges define the membrane of the region. Links between each membrane are used to pass Myxobacterial P systems between regions. This system is evolved by the execution of multi-environment P systems as outlined in Table 1. Myxobacteria in region \(i + 1, j + 1\) have moved into region \(i,j + 1\) using slime detection. The population of Myxobacteria in region \(i, j - 1\) spontaneously reversed direction from facing south to facing north. Population \(i - 1, j\) attempts to invade the populated region \(i, j\). By doing so, it causes a head–head collision, thus activating the C-signal pathway, causing a reversal in both P systems. During this time, the population of \(i + 1, j\) moved in besides the population of \(i, j\).

Fig. 1
figure 1

An illustrative representation of a multi-environment P system. Each square is a membrane region. Each tuple represents a Myxobacterial population P system

Full size image

3.2 Simulation examples

We present simulations of a Myxobacterial population as a colony. Default stochastic values are presented in Table 2 and draft Java code can be downloaded from https://github.com/acnash/Psystems. Video files of the simulations are available as online supplementary material. Each Myxobacterial P system was populated from a normal distribution with at most 400. In each region, the Myxobacteria were randomly assigned an initial direction.

Table 2 Default values for the stochastic selection of rules. Stochastic values are retrieved for valid rules (cellular conditions are correct), arranged descending by size and transformed onto the range (0, 1) for random selection
Full size table

3.2.1 Simulation examples of a Myxobacterial colony

A colony of up to 40,000 Myxobacteria was constructed about a \(40 \times 40\) multi-component P system substrate. The time step was a single evaluation of rules across the system and the simulation came to an end after 250 time steps had been evaluated. Three simulations were performed: first, a colony of Myxobacteria traversing the substrate under the influence of social and adventurous motility; second, a colony of Myxobacteria acting under the influence of adventurous behaviour; third, acting under the influence of social coordinated motility. The quantity of Myxobacteria at each multi-component P system region is represented by the integer values as specified in the Myxobacterial P system definition and is reflected by the peak high in Figs. 2, 3 and 4.

Fig. 2
figure 2

Simulation trajectory snapshots at 50 time intervals. S-motility and A-motility are enabled

Full size image
Fig. 3
figure 3

Simulation trajectory snapshots at 50 time intervals. The stochastic values for S-motility have seen set to 0, to ensure bacterial aggregation under complete adventurous motility

Full size image
Fig. 4
figure 4

Simulation trajectory snapshots at 50 time intervals. The stochastic values for A-motility have seen set to 0, to ensure bacterial aggregation under complete social motility

Full size image

During the early simulation steps of a Myxobacteria colony capable of social and adventurous exploration (Fig. 2), the Myxobacterial P system populations at the colony periphery orientated out towards the empty adjacent multi-component P system regions. Those Myxobacterial P systems within the heart of the colony would have potentially undergone a number of C-signal reversals from head-to-head collisions with immediate neighbours (\(r_{7}\)), spontaneous reversals (\(r_{10}\)) or exchange of Myxobacterial sub-populations between neighbouring Myxobacterial P systems (\(r_{{ 6 {\text{a}}}}\), \(r_{{ 6 {\text{b}}}}\) and \(r_{{ 6 {\text{c}}}}\)). Driven by the adventurous motility (\(r_{2}\)) of the periphery Myxobacterial P systems, immediate neighbours were able to detect local deposits of slime encouraging movement away from the colony. As the initial colony structure reduces, some of those adventurous and free-exploring Myxobacteria switched to S-motility (\(r_{9}\)), enabling neighbouring S-motility cells to begin mixing and moving amongst themselves (\(r_{{ 6 {\text{a}}}}\), \(r_{{ 6 {\text{b}}}}\) and \(r_{{ 6 {\text{c}}}}\)). Before too long, the periphery Myxobacterial P systems having broken away from the main colony to setup smaller colonies can switch back into A-motility (\(r_{8}\)) and once again venture out to form additional colonies.

By being able to turn towards the direction of slime deposits (\(r_{6a}\), \(r_{{ 6 {\text{b}}}}\) and \(r_{{ 6 {\text{c}}}}\)), switching direction randomly (\(r_{10}\)), and interacting with neighbouring adventurous P systems (\(r_{5}\)), there is a high degree of variability compared to the initial simulation parameters and an increased ability to form further sub-colonies. It is interesting to note, by decreasing the probability of a population switching from A-motility to S-motility, the colony immediately breaks up and begins randomly exploring the environment (data not shown). This is seen when the genes responsible for the Myxobacteria’s S-motility are knocked out, preventing the Myxobacteria from gripping onto neighbouring Myxobacteria [29].

A simulation of a Myxobacterial colony using only the A-motility (rules \(r_{2}\), \(r_{{ 3 {\text{a}}}} ,\; r_{{ 3 {\text{b}}}} ,\; r_{{ 3 {\text{c}}}}\), \(r_{{ 4 {\text{a}}}} ,\; r_{{ 4 {\text{b}}}} , \;r_{{ 4 {\text{c}}}}\)) rules for movement and removing the ability for the P systems to switch motility was performed for 250 time steps (Fig. 3). Compared to the normal-motility simulation, during the initial simulation steps the centre of the colony remains intact whilst showing only minor indication of the periphery Myxobacterial P system orientating towards the empty substrate to forward (\(r_{2}\)) into the empty substrate. Without S-motility, the Myxobacterial P systems are unable to maintain a bacterial aggregate, rather they move out into explore the empty substrate (\(r_{2}\)) or follow the slime deposits of Myxobacterial P system directly in front (\(r_{{ 3 {\text{a}}}}\)), or at an angle (\(r_{{ 3 {\text{b}}}}\), \(r_{{ 3 {\text{c}}}}\)). The eight orientations available can be clearly seen from step 100 to step 250, with a less dense covering of Myxobacteria as populations detected neighbouring slime on the angle off from their long axis.

Similarly, as with A-motility, when the movement of the population has been reduced to only S-motility (\(r_{1}\), \(r_{{ 6 {\text{a}}}}\), \(r_{{ 6 {\text{b}}}}\), and \(r_{{ 6 {\text{c}}}}\)) (Fig. 4), the eight possible orientations can also be discerned across the Myxobacterial P system population. Without adventurous motility, aggregates of P system Myxobacteria migrate along these orientations with little in the way of break-away Myxobacteria, where each colony having formed due to collisions (\(r_{7}\)) and spontaneous reversal (\(r_{10}\)).

To discern spatial differences over time between the three means of motility, the Myxobacterial P System distribution was compared to a uniform Poisson point null model, often called the complete spatial randomness (CRS) model [30]. At each time step, a quadrat count over the multi-environment P system space was calculated. As P system distributions become more random, the quadrat count calculation moves closer to accept the null. Calculations were performed per step over the three motilities and a log transformation of their p value was visualised (Fig. 5). The initial colony structure results in an identical transformed p value over the three models. As each model evolves the transformed p value increases towards the null model (p > 0.05). The adventurous motility, as per the visible ordered distribution in Fig. 3, retains the lowest p value. As the simulations evolve, social motility shows the greatest degree of randomness (as per Fig. 4). Combined, the transformed p value of normal motility is neither as random as just social motility nor as ordered as adventurous motility.

Fig. 5
figure 5

Log transformation of p values of quadrant counting tests calculations per step

Full size image

4 Conclusion

Using a multi-component P system substrate and a model P system to represent dynamic quantities of Myxobacteria, this work presented simulations of a swarming and aggregating Myxobacterial population. The results demonstrate how both modes of bacterial motility, social and adventurous were necessary to replicate both swarming of bacteria from the colony edge, exploration of the empty substrate and intermittent sub-colony formation. A simulation using only social behaviour explored the substrate by forming uniformed sub-colonies, yet with no means of acting as solitary agents, the distribution of populations was along the eight orientations. Whereas, adventurous motility only simulations results in Myxobacterial P systems at the colony periphery venturing into the empty substrate and in turn leads to neighbouring P systems breaking from the colony. Yet, with no mixing of P system populations (from social behaviour), had the simulation been left indefinitely, it is quite possible that the complete population would have been distributed equally over the substrate space.

With all languages being Turing equivalent, it is possible to implement a P System model using any language as mathematically, all programming languages are Turing equivalent. However, Java was selected given its native support for encapsulation and inheritance. Object-orientated languages such as Java and C++ lend themselves very well to the paradigm of membrane computing. Although of interest are the efforts to encapsulate P system rules in a formal language [31], and very recently using in vitro techniques to produce a real P system [32].

Despite the swarming and aggregating dynamics observed using only a small set of P system rewrite rules, there are a number of limitations; orientation is limited to the eight adjacent spaces in the multi-environment P system substrate; the time step is a simple count and the execution of rules are determined stochastically with no consideration of temporal dynamics; and, slime remains indefinitely, reducing the A-motility to a random walk as the simulation progresses. Compared with recent particle continuum [33] and continuous [34] based modelling, conventional models of Myxobacteria continue to dominate, chiefly through their ability to model structural properties of bacterium. Finally, there is a degree of abstraction from the traditional cell-line P system paradigm in the way in which the mixing of bacteria is treated. Traditionally, rewrite rules and the diffusion of cell content would be required to combine Myxobacterial populations. However, given the framework presented in this study, further exploration into locally induced global phenomena of bacterial colonies using P systems is entirely possible.