Spatial biology of Ising-like synthetic genetic networks

Background Understanding how spatial patterns of gene expression emerge from the interaction of individual gene networks is a fundamental challenge in biology. Developing a synthetic experimental system with a common theoretical framework that captures the emergence of short- and long-range spatial correlations (and anti-correlations) from interacting gene networks could serve to uncover generic scaling properties of these ubiquitous phenomena. Results Here, we combine synthetic biology, statistical mechanics models, and computational simulations to study the spatial behavior of synthetic gene networks (SGNs) in Escherichia coli quasi-2D colonies growing on hard agar surfaces. Guided by the combined mechanisms of the contact process lattice simulation and two-dimensional Ising model (CPIM), we describe the spatial behavior of bi-stable and chemically coupled SGNs that self-organize into patterns of long-range correlations with power-law scaling or short-range anti-correlations. These patterns, resembling ferromagnetic and anti-ferromagnetic configurations of the Ising model near critical points, maintain their scaling properties upon changes in growth rate and cell shape. Conclusions Our findings shed light on the spatial biology of coupled and bistable gene networks in growing cell populations. This emergent spatial behavior could provide insights into the study and engineering of self-organizing gene patterns in eukaryotic tissues and bacterial consortia. Supplementary Information The online version contains supplementary material available at 10.1186/s12915-023-01681-4.


Supplementary Note 1: Simulation of the Ising model during the growth of a quasi 2D colony
To study the effects of Ising-like gene-to-gene interactions in a population of bacterial cells growing as a quasi-two-dimensional (2D) colony, we created a lattice model that combines spatial growth dynamics and the Ising model from statistical mechanics to represent the cell state differentiation induced by our genetic network.We combined the simulation of the two-dimensional Ising Model (IM), using the Metropolis algorithm [37], with the Contact Process (CP) lattice model of spatial stochastic growth [36] to represent cell population dynamics of the bacterial colonies observed in our range expansion experiments.We name this model as CPIM (Contact Process -Ising Model).CPIM consists of an interacting particle system that follows colonization, extinction process mixed with differentiation dynamics of particles on a two-dimensional lattice (L) of N sites.An important consideration is that colonization/extinction (or birth/death) processes operate at the scale of the entire bacterial cell while differentiation (induction of gene expression) occurs at the scale of molecular networks within living cells.We use this separation of organizational scales to simplify the mixing of the CP and IM.
Each site of this lattice, representing square sectors of the agar surface, can be in one of four states S = {∅, * , −1, +1}, which represent vacant locations (state ∅, white squares), locations occupied with undifferentiated cells (state * , black squares), and locations occupied with differentiated cells expressing RFP (state +1, magenta squares) or GFP (state −1, green squares) reporter genes respectively.The state of a site x at a time t can be described as ξ t (x) [38], an interacting particle system ξ t : L → S that, for each time point t, assigns one state from S to all sites in the lattice.This spatial stochastic process is determined by the following lattice reactions from three different sets of processes: i) Birth, death, and differentiation processes: cells die (lattice squares suffer local extinction) and divide at rates δ and b, respectively.For this process an event (local extinction) occurring at rate δ, means that the times ∆t i between successive occurrences (lifespan of the occupied states) have an exponential distribution with parameter δ; that is, P(∆t i ≤ 1) = 1 − exp(−δt).Colonization reactions (utilization of vacancy, state ∅) all depend on the state of nearby sites as they are the sources of colonization propagules.Grouping these states as a super-state Ω = { * , −1, +1} we have that ξ t : L → {∅, Ω} corresponds to a CP.Vacant sites can be colonized only by Nearest Neighbors (NN; ∥x − y∥ ≤ 1 ∀x, y ∈ L; i.e. north, south, east and west), originating an advancing (2-D) front of growing cells that colonize available habitat in the landscape.This spatial spread, for all ω ∈ Ω, is represented by The rate of colonization (bρ ω ) is proportional to both, the birth rate b and the local density ρ ω of ω ∈ Ω states in the NN neighborhood of the focal vacant (∅ x ) site x ∈ L.
ii) Differentiation and gene expression process.In our experimental system, there is always the possibility that cells do not express our synthetic gene network.We model this as an undifferentiated state * and therefore at the beginning of most simulations, a square representing an undifferentiated cell ( * state) is placed in the center of the lattice.Then, we proceed by a Monte Carlo method repeating as many times as there are sites on the lattice to produce one generation (step) of our simulation: a site on the lattice is randomly selected, and if this focal site is vacant (∅ state), a random neighbor is selected.If this neighbor is occupied (Ω states, occurring with probability, ω∈Ω ρ ω ), it colonizes the vacant site with probability bρ ω • dt.On the other hand, if the focal site is occupied by an undifferentiated state ( * ), if it survives (with probability 1 − δ • dt), it can differentiate into an RFP-expressing (+1) or a GFP-expressing (−1) state with a certain probability α, which is represented by * As members of the CP's occupancy class Ω, these differentiated states can also die or give birth with rates δ and bρ * respectively.Notice that, unlike the birth-death process, internal cell states changes, like genetic differentiation, are conditional upon cell survival, and instead of being determined by rates, they are ruled by probabilities which are conditional on a particle's survival.We use this separation of organizational scales between a complete bacterial cell and a small subset of its constituent components (our synthetic gene network) to make the CP model of vacancy utilization independent from our different distinctions of occupancy (Ω states).
iii) Ising-like cellular states alignment process: as our genetic network consists of a bi-stable switch locally coupled by diffusing small molecules (auto-inducers) that trigger flips between two reporter gene states (GFP and RFP); we model it as spin variables and represent them using the IM (states -1 and +1).To determine if lattice squares representing differentiated cells can change from one state to another (RFP to GFP or GFP to RFP) we use the Metropolis algorithm.Thus, when a site in the lattice that is occupied by a differentiated state (-1 or +1) is randomly selected during a Monte Carlo step, if the square survives extinction (with probability 1 − δ • dt) the following computation takes place: We calculate the interaction energy of the focal site (i = x ∈ L) and its neighbors (j = y ∈ N ) using H = −J ⟨ij⟩ σ i σ j (Eq.1; main text).Here, σ i and σ j ∈ {−1, +1} represent spin variables and the notation ⟨ij⟩ indicates that the sum is between neighboring pairs of sites, allowing only a short-range interaction.In most simulations, we considered neighborhoods N to be NN but we also explored the Next Nearest Neighbors (NNN) case (∥x − y∥ ≤ 2 ∀x, y ∈ L; where ∥ • ∥ is the Manhattan metric).Next, we determine the change in the interaction energy ∆H, which results if the spin state at that site is flipped.For the interaction energy calculations, we assumed that J = +k B (where k B is the Boltzmann constant) for ferromagnetic systems and J = −1k B for anti-ferromagnetic systems (see genetic constructs; main text).This assumption implies that the interaction energy between neighboring sites is minimized when they have the same state in a ferromagnetic system, or the opposite state in an anti-ferromagnetic system (Eq.1; main text).Thus, if ∆H < 0, the change in the spin state ([−/+] x P − → [+/−] x ) is always accepted (probability P = 1) since this change decreases the total energy of the system.If ∆H ≥ 0, which means an increase in the total energy of the system, the change in the spin state is accepted only with probability P = e −∆H/k B T [39].In the classical IM, T is the absolute temperature and corresponds to the control parameter of the system.This means that varying T determines the transition between ordered and disordered configurations of the system.In our model, T represents a parameter that determines the coupling strength between cells in the experimental system: a small value of T represents a strong coupling, while a large value represents a weak coupling.

Simulations of the model
In this model, one generation time corresponds to the moment when statistically all the sites in the lattice have had the possibility of being visited.In a 256x256 lattice, we run simulations for 6800 to 9800 generation times.Without loss of generality, we consider k B = 1.

Considerations of our Model assumptions
It is also important to note that a lattice site in the CPIM simulations can have two interpretations.As described above, a lattice site in the CPIM simulations can represent a cell.However, a lattice site in this model can also represent a cell territory.Considering lattice sites as cell territories is particularly relevant since it allows to map a lattice site to a sector in a colony.In this interpretation, besides considering lattice sites as cell territories (and not individual cells), the death rate represents the clearance rate of territories, and the division/birth rate represents the local colonization rate.However, the usefulness of the model remains the same: local interactions can give rise to large-scale order differences.

Temporal considerations
The CP has rates defined by the Poisson distribution.In the Mean Field (MF), the proportion of occupied sites, p, can be modeled as, ṗ = bp(1−p)−δp and at equilibrium we have p = 1 − δ/b.Thus, b > δ is needed for the system to persist, having a residual vacancy which corresponds to the time units, b −1 , it takes particles to replicate, times their dead/extinction rate.In the case of our colonies, local extinction (death) is negligible at the scale of our observations and compared with the colony expansion rate.Thus in most simulations, we used dt ≡ 1 with b = 0.03 and δ = 0.00001.This means, on average, it takes 100 generations in our simulations for 3 particles to replicate when provided with an adjacent vacancy and 100000 generations for them to cease to exist.This mostly takes place at the expansion front of the simulated colonies.The CP representing our colonies is very dense (≈ 99, 9% occupancy, Additional file 3, Fig. S13A) as b/δ >> 1.We decide to use the maximum possible rate for the spin flips of the IM at the timescale of one generation (dt).
On infinite temperature, T → ∞, all spins flip on average in one generation (Monte Carlo step; dt ≡ 1).On the contrary, at zero temperature, spins flip with probability one if the energy difference is negative and not at all if it is positive.When the temperature is low, spins always flip for negative energy difference and with exponentially decaying probabilities when energy differences are positive.
The important quantity when simulating the IM is the energy difference obtained by flipping a state ξ t (x) ∈ {−1, +1} to its opposite sign.Here N + and N − correspond to the number of sites with spin states +1 and −1 in the neighborhood N of x respectively.The Metropolis algorithm accepts a spin flip with probability λ flip depending on the energy difference between states before and after the proposed flip (see Additional file 3, Fig. S17).This probability is not normalized by the partition function and all configurations with ∆H ≤ 0 are flipped within a Monte Carlo step (λ flip = 1).On the contrary, if ∆H > 0 then we flip with probability (1) Notice that the Markov process produced by a Monte Carlo step on the lattice (see Additional file 1, Fig. S1C) is not meant to represent flipping rates but rather to draw configurations that are consistent with the Boltzmann distribution given by the Hamiltonian of Eq.1 (main text).Therefore in the CPIM, energy-favorable flips are considered to be an instantaneous process that is conditional to CP's particle survival.Without loss of generality, as it only introduces transient states, we include an undifferentiated state which turns into spins (differentiated states) with no bias with probability α = 0.1.Thus, in one generation of the CPIM, 10% of surviving undifferentiated occupied sites, become spin variables of the IM.

Majority rules
The CPIM has transition probabilities that depend on the configurations of neighborhoods.The colonization of an empty site in the CP occurs at rate b when the vacant site is fully surrounded by occupied ones (or probability b • dt).In Additional file 3, Fig. S17A we see such a vacant site fully surrounded (ρ * = 1) by the undifferentiated state * .On the contrary, sites can go extinct independently of neighborhood configuration at rate δ (or probability δ • dt).When making distinctions of occupancy states (Ω), their neighborhood densities ρ ω∈Ω determine their colonization likelihood.For the IM, there are a finite number of neighborhood configurations of sites that are in spin states.As we see in Additional file 3, Fig. S17B and S17C and Equation 1 there is an even smaller number of possible spin occupancy patterns (N − , N + ) for these neighborhoods.For the case of NN interactions, if Temperature is T = 1 and J = k −1 B , we only have 11 possible values of ∆H and their corresponding flipping probabilities: That is, even at small temperature values, all spin-flips with favorable energy are executed within the time scale of one generation on all surviving particles.As we do not have detailed information on the diffusion of HSL molecules in the agar, nor on their transport rates we decide to use the IM as implemented by the Metropolis algorithm when hybridizing it with the CP.Most likely configurations with higher energy differences in favor of a flip would happen at a higher rate, the large order of magnitudes in differences will introduce numerical complexity when running simulations.Thus, as a first approach, we decide to leave the IM as it is implemented in the Metropolis algorithm and use the flipping pseudo rates as rates of our model.In practice, this means that a few neighboring sites favoring a given flip will be enough to trigger it.If vacancy is abundant, however, the CP alters the IM as persistent isolated sites (surrounded by vacancy) have zero interaction energy regardless of their spin value and will flip with probability one.Frustrated sites where the spin occupancy N + − N − = 0 also flip instantaneously.This reflects the bi-stable nature of the genetic network and our ignorance of its characteristic time scales (see Section 0.3 hereunder) 0.2 Supplementary Note 2: Finite-size scaling analysis Studying the effects of the finite size of the system on observables that characterize the model's Ising-like critical behavior such as the average absolute value of the magnetization per site ⟨|M |⟩ (Additional file 3, Fig. S13B), and consequently, the magnetic susceptibility per site χ = N (⟨|M | 2 ⟩ − ⟨|M |⟩ 2 ) (Additional file 3, Fig. S13C), allows us to assess whether the phase transition exhibited in the CPIM falls under the universality class of that of the Ising model.This can be done by determining the critical exponents ν, γ, and β (Additional file 3, Fig. S13E, S13F, and S13G respectively; Additional file 4, Table S6) as done in Weber M, Buceta J., 2016) [22].This task was carried out by running CPIM simulations for a range of temperatures T ∈ [1.0, 3.6], on different lattice sizes L = {16, 32, 64, 128}, using the same model parameters b = 0.03, d = 0.00001, and α = 0.1.All simulations were initialized from a single undifferentiated lattice site, run for 50, 000 generations to ensure the stabilization of the colonization/extinction dynamics (Additional file 3, Fig. S13A), and finally, data was recorded for 15, 000 generations.Each simulation was repeated 10 times (replicates), critical exponents are means over the replicates and their respective standard deviations are shown.
As seen in Additional file 3, Fig. S13B and S13C, the critical temperature value T c at which the (magnetic) phase transition occurs is difficult to pinpoint when the system's size is too small.This problem can be circumvented by means of the cumulant intersection method developed by K. Binder [52], which consists of comparing the Binder cumulant B (Equation 2) as a function of the temperature T for different lattice sizes L (Additional file 3, Fig. S13D).
The intersection of the Binder cumulant B curves marks the critical temperature of the system in the thermodynamic limit T c ≈ 2.27 ± 0.01 (Additional file 3, Fig. S13D, dashed line) since B does not depend on the system size at the critical point [53].
After finding the critical temperature T c at which the system undergoes an Isinglike phase transition, a fit to the finite-size scaling relation was used to determine the value of the critical exponent ν: ν = 2.04 ± 0.02, where T N c corresponds to the position of the peak of the magnetic susceptibility per site for each lattice size studied (Additional file 3, Fig. S13E) and f is a constant.
Furthermore, both the height of the peak of the magnetic susceptibility per site χ N max and the average magnetization per site at the critical point ⟨|M |⟩ N c depend on the system's size following power-laws (Additional file 3, Fig. S13F and S13G respectively), so a fit to where g is a constant, yields γ = 1.79±0.05.Finally, fitting the simulations data to where h is a constant, reveals the critical exponent β: β = 0.16 ± 0.02.All critical exponents calculated in this work are means (of generation time averages) over 10 replicates, and standard deviations of these means are presented as uncertainty measures.The critical exponents calculated here and their respective values for the Weber-Buceta model [22] and a numerical solution of the two-dimensional Ising model are summarized in Additional file 4, Table S6).

Supplementary Note 3: Time series analysis
To gain insight into the dynamics of reporter gene expression patterns in our colonies and how to relate them with the spin flipping rates produced by the Metropolis algorithm of the Ising Model, we compared pictures of experimental colonies taken at three snapshots separated by 4 hours of growth.Examples of these colonies are shown in Additional file 3, Fig. S15.These correspond to: (i) 14 hours, (ii) 18 hours, and (iii) 22 hours post-inoculation respectively.In Additional file 3, Fig. S18 we show how we can estimate colony expansion by measuring changes in its area and diameter.Notice that colony expansion is significantly slower when comparing 18 hours and 22 hours, in contrast to what is observed when comparing 14 hours with 18 hours (see Additional file 3, Fig. S16A and S16B).By comparing colonies between two successive snapshots, we can define a core and an expansion ring as shown in Additional file 3, Fig. S16C.We can estimate the sectors of the colony which change the pattern of expression and observe, as shown in Additional file 3, Fig. S16D, that the core, defined by two consecutive snapshots, reflects more similarity in the pattern of gene expression (white pixels) as the colonies grow older.This can be due to entry into stationary phase or other stresses experienced by the colony as its radius becomes larger.
Although our CPIM approach considers a constant and spatially homogeneous rate of spin flipping, an interesting avenue for future theory development is to develop models where this can vary in space and time.To be able to develop such models, it would be extremely important to have experiments that measure the dynamics of these patterns in gene expression under different ecological conditions of food supply and waste removal.

Supplementary Note 4: Comparisons with the work of Weber and Buceta
The work of Weber and Buceta (WB; [22]) is a relevant study to consider when analyzing our experimental results.In their work, they propose a proof of concept model of a gene regulatory circuit.Although our theoretical analysis and simulations differ considerably from their work, both are inspired by the Ising model and can be used in a complementary fashion to interpret our experiments.The first consideration to consider is that our work is guided by experiments and we develop theory to interpret those specific experiments.Secondly, our experiments are on range expansion of an E. coli colony growing on a two-dimensional hard agar surface.Therefore, we focus on a spatially structured model of growth and not on modeling the specifics of the underlying molecular network.However, we do use a complementary model to calibrate the molecular constructs we use for coupling our genetic switch (see Eq. 2 and Fig. 2.; main text).Although in both cases the use of autoinducer molecules is suggested to couple bi-stable switches, our coupling constructs are different from the proof of concept proposed in WB.Another important difference is that in our work we explored both, ferromagnetic and anti-ferromagnetic interactions experimentally.
Although our simulations, as well as genetic details of our molecular constructs, and those of WB are different, we can interpret our experiments while comparing both theoretical approaches (see Discussion; main text).The switch studied in [23] corresponds to proteins LacI (U in WB) and λ (V in WB) which define a concentration pair (u, v) that acts as a toggle switch which in ideal theoretical conditions can be thought to exists in two alternative configurations (phenotypes) s 1 (u s1 , v s1 ) = (5.38,0.508)nM (spin up, state +1) or s 2 = (u s2 , v s2 ) = (v s1 , u s1 )nM (spin down, state -1).As they point out, in a volume of 1.5µm 3 (typical E. coli cell) this corresponds to 0.5 and 5 molecules per cell on average.A number that will be challenging to measure experimentally.Even considering noisy fluctuations with an amplitude reaching 15 molecules per cell (Fig. 2 in WB), it would still be a challenging experiment.In our case, we used fluorescent reporters which are expressed in concentrations easily observable by fluorescent microscopy, while focusing our attention on a macroscopic colony expanding over a surface of hard agar.WB, on the other hand, examine a different demographic context when they turn their attention to an unstructured theoretical population (Fig. 5 therein, [22]).Spatial structure is of crucial importance in our experimental setup and therefore we use the 2D Contact Process (CP) to model our expanding colony (see main text, and Sections 0.1 and 0.3 of this document).
Our experimental genetic switch is also different from the switch proposed in the proof of concept circuit in WB, which is based on the toggle switch proposed by [23].In their simulations, WB uses an original toggle switch (module 1 in their Figure 2), in which repressors (LacI and λcI) are regulated by constitutive-repressible promoters, in parallel to a HSL-regulated expression of receptors (LuxR/LasR) and a second copy of repressors (LacI and λcI) (module 2 and 3 in their Figure 2).Therefore, they have HSL-regulated receptors (LuxR/LasR) and two copies of repressors.In our system, we have only one copy of each repressor.Also, the LuxR/LasR are expressed constitutively.The coupling signals in our system regulate the expression of inducerproducing enzymes LuxI and LasI and repressors LacI and TetR; whereas in WB, HSL regulates one of the two copies of repressors and LuxR/LasR.Thus, we have inducible-repressible promoters that are induced by HSL and repressed by repressors.These promoters regulate the repressors, LuxI/LasI, and the reporters.In our system, repressors (LacI/TetR) and inducer (LuxI/LasI) genes are regulated at the same time by these inducible-repressible promoters (pLux-pTet regulating LacI and LuxI, pLas-pLac regulating TetR and LasI).We use LacI and TetR as repressors, while they use LacI and λcI.
WB assume no cross-talk while we do observe it to exist experimentally at high auto-inducer concentrations.This fact, together with the induction of stationary phase, limits the time of growth and the size of colonies to be considered in our experiment.Partition of transcription space in E.coli means that the polymerase changes its affinity for some promoters at different phases.
WB develop an "effective" approach by focusing on the fact that U and V act as their own auto-inducers.They use this fact to formulate a model which depends on an effective transport parameter D regulating concentrations of actors in intracellular (U i and V i ) as well as extracellular compartments (U e and V e ).Transport between compartments is modeled as rD where r represents the ratio between compartment volumes, r = V cell /V ext .A crucial point is their assumption that diffusion across compartments is at least one order of magnitude different than diffusion of autoinducers in the extracellular compartment.The CPIM does not attempt to model transport or concentrations of proteins as it rather focuses on spatial occupancy.It models the spatial spread of occupancy over vacancy.To model the genetic switch we simply make three distinctions of occupancy types (lattice sites): un-differentiated site (state * ), spin down (state -1), and spin up (state +1).We propose our model to be also an "effective" theory in the sense that our experimental spatial data (pixels of a fluorescent image) corresponds to area sectors of agar, where autoinducer molecules that have accumulated in both intracellular as well as extracellular space can induce expression of our reporters in that territory.In the case of ferromagnetic experiments a pixel induces the same state on their neighbors and in the case of anti-ferromagnetic case, it induces the opposite state.We represent this by the Ising Hamiltonian and its corresponding sign for parameter J values.Only our ferromagnetic experiments can be interpreted in the light of WB.However, we do not have detailed information about how effective are autoinducers when triggering changes in gene expression, nor their actual concentrations in the agar surface or their diffusion rates.Our model is as generic as possible to be able to be calibrated when further experiments in microfluidic devices are performed.
Another point to consider is that when the reproductive number of the Contact Process is high enough, our model exhibits exponents that are more similar to the IM than those of WB.As in the IM, in CPIM, temperature T and coupling J are independent.WB assume a concentrated culture condition to develop an order parameter.We use the temperature T parameter of the IM as our control parameter.The most important part is the analogy between temperature and auto-inducer concentration and transport.This produces a critical value of inducer which triggers the transition.Biological parameters that characterize the binding of the inducers C6HSL and C12HSL to the RFP and GFP promoters in cells carrying the reporter vector (Table S1), the ferromagnetic system (Table S2) or the anti-ferromagnetic system (Table S3).The values were obtained from the fit to Eq. 2 of the data of red and green fluorescent protein synthesis rates of liquid cultures of cells carrying these systems as a function of the inducer concentrations.Power-law exponents of ferromagnetic populations simulated with CPIM at T = Tc, ferromagnetic colonies of rod-shaped and spherical cells (Fig. 4d), and ferromagnetic colonies of spherical cells grown in minimal medium supplemented with glucose (GLU) or glycine (GLY) (Fig. 4e).The exponents were estimated using the r plfit(k,'hist') algorithm [69] or by the least squares method.

Table S4 :
Scaling exponent γ of populations with ferromagnetic interactions.

Table S5 :
Scaling exponent γ of simulated populations with ferromagnetic interactions at different birth rates.

Table S6 :
Numerical critical exponents and critical temperature of the two-dimensional Ising model, Weber-Buceta model, and the Contact Process Ising Model (CPIM).Data of the Weber-Buceta model and the 2-D Ising model was taken from Weber M, Buceta J.,2016 [22]and E. Ibarra-García-Padilla et al., 2016 [40]respectively.

Table S7 :
Plasmids used in this study