View source for BIO-LGCA

In [[Computational biology|computational]] and [[Mathematical and theoretical biology|mathematical biology]], a '''biological lattice-gas cellular automaton''' (BIO-LGCA) is a discrete model for moving and interacting biological agents,<ref>{{Cite journal|last1=Deutsch|first1=Andreas|last2=Nava-Sedeño|first2=Josué Manik|last3=Syga|first3=Simon|last4=Hatzikirou|first4=Haralampos|date=2021-06-15|title=BIO-LGCA: A cellular automaton modelling class for analysing collective cell migration|journal=PLOS Computational Biology|language=en|volume=17|issue=6|pages=e1009066|doi=10.1371/journal.pcbi.1009066|issn=1553-7358|pmc=8232544|pmid=34129639|bibcode=2021PLSCB..17E9066D |doi-access=free }}</ref> a type of [[cellular automaton]]. The BIO-LGCA is based on the [[Lattice gas automaton|lattice-gas cellular automaton]] (LGCA) model used in fluid dynamics. A BIO-LGCA model describes cells and other motile biological agents as point particles moving on a discrete lattice, thereby interacting with nearby particles. Contrary to classic cellular automaton models, particles in BIO-LGCA are defined by their position and velocity. This allows to model and analyze active fluids and collective migration mediated primarily through changes in momentum, rather than density. BIO-LGCA applications include cancer invasion<ref>{{Cite journal|last1=Reher|first1=David|last2=Klink|first2=Barbara|last3=Deutsch|first3=Andreas|last4=Voss-Böhme|first4=Anja|date=2017-08-11|title=Cell adhesion heterogeneity reinforces tumour cell dissemination: novel insights from a mathematical model|journal=Biology Direct|volume=12|issue=1|pages=18|doi=10.1186/s13062-017-0188-z|issn=1745-6150|pmc=5553611|pmid=28800767 |doi-access=free }}</ref> and [[Cancer|cancer progression]].<ref name=":0">{{Cite journal|last1=Böttger|first1=Katrin|last2=Hatzikirou|first2=Haralambos|last3=Voss-Böhme|first3=Anja|last4=Cavalcanti-Adam|first4=Elisabetta Ada|last5=Herrero|first5=Miguel A.|last6=Deutsch|first6=Andreas|date=2015-09-03|editor-last=Alber|editor-first=Mark S|title=An Emerging Allee Effect Is Critical for Tumor Initiation and Persistence|journal=PLOS Computational Biology|language=en|volume=11|issue=9|pages=e1004366|doi=10.1371/journal.pcbi.1004366|issn=1553-7358|pmc=4559422|pmid=26335202|bibcode=2015PLSCB..11E4366B |doi-access=free }}</ref>

== Model definition ==
As are all cellular automaton models, a BIO-LGCA model is defined by a lattice <math>\mathcal{L}</math>, a state space <math>\mathcal{E}</math>, a neighborhood <math>\mathcal{N}</math>, and a rule <math>\mathcal{R}</math>.<ref name=":1">{{Citation|title=Mathematical Modeling of Biological Pattern Formation|url=http://dx.doi.org/10.1007/0-8176-4415-6_3|work=Cellular Automaton Modeling of Biological Pattern Formation|series=Modeling and Simulation in Science, Engineering and Technology|year=2005|pages=45–56|place=Boston, MA|publisher=Birkhäuser Boston|doi=10.1007/0-8176-4415-6_3|isbn=978-0-8176-4281-5|access-date=2021-05-25}}</ref>

* The lattice (<math>\mathcal{L}</math>) defines the set of all possible particle positions. Particles are restricted to occupy only certain positions, typically resulting from a regular and periodic [[Tessellation|tesselation]] of space. Mathematically, <math>\mathcal{L}\subset\mathbb{R}^d</math> is a discrete subset of <math>d</math>-dimensional space.
* The state space (<math>\mathcal{E}</math>) describes the possible states of particles within every lattice site <math>\mathbf{r}\in\mathcal{L}</math>. In BIO-LGCA, multiple particles with different velocities may occupy a single lattice site, as opposed to classic cellular automaton models, where typically only a single cell can reside in every lattice node simultaneously. This makes the state space slightly more complex than that of classic cellular automaton models (see below).
* The neighborhood (<math>\mathcal{N}</math>) indicates the subset of lattice sites which determines the dynamics of a given site in the lattice. Particles only interact with other particles within their neighborhood. Boundary conditions must be chosen for neighborhoods of lattice sites at the boundary of finite lattices. Neighborhoods and boundary conditions are identically defined as those for regular cellular automata (see [[Cellular automaton#Overview|Cellular automaton]]).
* The rule (<math>\mathcal{R}</math>) dictates how particles move, proliferate, or die with time. As every cellular automaton, BIO-LGCA evolves in discrete time steps. In order to simulate the [[system dynamics]], the rule is synchronously applied to every lattice site at every time step. Rule application changes the original state of a lattice site to a new state. The rule depends on the states of lattice sites in the interaction neighborhood of the lattice site to be updated. In BIO-LGCA, the rule is divided into two steps, a probabilistic interaction step followed by a deterministic transport step. The interaction step simulates reorientation, birth, and death processes, and is defined specifically for the modeled process. The transport step translocates particles to neighboring lattice nodes in the direction of their velocities. See below for details.

=== State space ===
[[File:Hexnode.png|thumb|The substructure of a BIO-LGCA lattice site with six velocity channels (corresponding to a 2D hexagonal lattice) and a single rest channel. In this case <math>b=6</math>, <math>a=1</math>, and the carrying capacity <math>K=7</math>. Channels 2, 3, 6 and 7 are occupied, thus the lattice configuration is <math>\mathbf{s}=(0,1,1,0,0,1,1)</math>, and the number of particles is <math>n\left(\mathbf{s}\right)=\sum_{i=1}^Ks_i=4</math>.]]
For modeling particle velocities explicitly, lattice sites are assumed to have a specific substructure. Each lattice site <math>\mathbf{r}\in\mathcal{L}</math> is connected to its neighboring lattice sites through vectors called "velocity channels", <math>\mathbf{c}_i</math>, <math>i\in\{1,2,\ldots,b\}</math>, where the number of velocity channels <math>b</math> is equal to the number of nearest neighbors, and thus depends on the lattice geometry (<math>b=2</math> for a one-dimensional lattice, <math>b=6</math> for a two-dimensional hexagonal lattice, and so on). In two dimensions, velocity channels are defined as <math>\mathbf{c}_i=\left(\cos\frac{2\pi i}{b},\sin\frac{2\pi i}{b}\right)</math>. Additionally, an arbitrary number <math>a</math> of so-called "rest channels" may be defined, such that <math>\mathbf{c}_i=(0,0)</math>, <math>i\in\{b+1,b+2,\ldots,b+a\}</math>. A channel is said to be occupied if there is a particle in the lattice site with a velocity equal to the velocity channel. The occupation of channel <math>\mathbf{c}_i</math> is indicated by the occupation number <math>s_i</math>. Typically, particles are assumed to obey an [[Pauli exclusion principle|exclusion principle]], such that no more than one particle may occupy a single velocity channel at a lattice site simultaneously. In this case, occupation numbers are Boolean variables, i.e. <math>s_i\in\mathcal{S}=\{0,1\}</math>, and thus, every site has a maximum [[carrying capacity]] <math>K=a+b</math>. Since the collection of all channel occupation numbers defines the number of particles and their velocities in each lattice site, the vector <math>\mathbf{s}=\left(s_1,s_2,\ldots,s_{K}\right)</math> describes the state of a lattice site, and the state space is given by <math>\mathcal{E}=\mathcal{S}^K</math>.

=== Rule and model dynamics ===
The states of every site in the lattice are updated synchronously in discrete time steps to simulate the model dynamics. The rule is divided into two steps. The probabilistic interaction step simulates particle interaction, while the deterministic transport step simulates particle movement.

==== Interaction step ====
Depending on the specific application, the interaction step may be composed of reaction and/or reorientation operators.

The reaction operator <math>\mathcal{A}</math> replaces the state of a node <math> \mathbf{s}</math> with a new state <math>\mathbf{s}^{\mathcal{A}}</math> following a [[Markov chain|transition probability]] <math>P\left(\left. \mathbf{s}\rightarrow \mathbf{s}^{\mathcal{A}}\right| \mathbf{s}_{\mathcal{N}} \right)</math>, which depends on the state of the neighboring lattice sites <math>\mathbf{s}_{\mathcal{N}}</math> to simulate the influence of neighboring particles on the reactive process. The reaction operator does not conserve particle number, thus allowing to simulate birth and death of individuals. The reaction operator's transition probability is usually defined ''ad hoc'' form phenomenological observations.

The reorientation operator <math>\mathcal{O}</math> also replaces a state <math>\mathbf{s}</math> with a new state <math>\mathbf{s}^{\mathcal{O}}</math> with probability <math>P\left(\left. \mathbf{s}\rightarrow \mathbf{s}^{\mathcal{O}}\right| \mathbf{s}_{\mathcal{N}} \right)</math>. However, this operator conserves particle number and therefore only models changes in particle velocity by redistributing particles among velocity channels. The transition probability for this operator can be determined from statistical observations (by using the [[Principle of maximum caliber|maximum caliber principle]]) or from known single-particle dynamics (using the discretized, steady-state angular probability distribution given by the [[Fokker–Planck equation|Fokker-Planck equation]] associated to a [[Langevin equation]] describing the reorientation dynamics),<ref>{{Cite journal|last1=Nava-Sedeño|first1=J. M.|last2=Hatzikirou|first2=H.|last3=Peruani|first3=F.|last4=Deutsch|first4=A.|date=2017-02-27|title=Extracting cellular automaton rules from physical Langevin equation models for single and collective cell migration|url=http://dx.doi.org/10.1007/s00285-017-1106-9|journal=Journal of Mathematical Biology|volume=75|issue=5|pages=1075–1100|doi=10.1007/s00285-017-1106-9|pmid=28243720|s2cid=32456636|issn=0303-6812}}</ref><ref>{{Cite journal|last1=Nava-Sedeño|first1=J. M.|last2=Hatzikirou|first2=H.|last3=Klages|first3=R.|last4=Deutsch|first4=A.|date=2017-12-05|title=Cellular automaton models for time-correlated random walks: derivation and analysis|url=http://dx.doi.org/10.1038/s41598-017-17317-x|journal=Scientific Reports|volume=7|issue=1|page=16952|doi=10.1038/s41598-017-17317-x|pmid=29209065|pmc=5717221|arxiv=1802.04201 |bibcode=2017NatSR...716952N |issn=2045-2322}}</ref> and typically takes the form <math display="block">P\left(\left. \mathbf{s}\rightarrow \mathbf{s}^{\mathcal{O}}\right| \mathbf{s}_{\mathcal{N}} \right)
=\frac{1}{Z}e^{-\beta H\left(\mathbf{s}_{\mathcal{N}}\right)}
\delta_{n\left(\mathbf{s}\right),n\left(\mathbf{s}^{\mathcal{O}}\right)}</math> where <math>Z</math> is a normalization constant (also known as the [[Partition function (mathematics)|partition function]]), <math>H\left(\mathbf{s}_{\mathcal{N}}\right)</math> is an energy-like function which particles will likely minimize when changing their direction of motion, <math>\beta</math> is a free parameter inversely proportional to the randomness of particle reorientation (analogous to the [[Thermodynamic beta|inverse temperature]] in thermodynamics), and <math>\delta_{n\left(\mathbf{s}\right),n\left(\mathbf{s}^{\mathcal{O}}\right)}</math> is a [[Kronecker delta]] which ensures that particle number before <math>n\left(\mathbf{s}\right)</math> and after reorientation <math>n\left(\mathbf{s}^{\mathcal{O}}\right)</math> is unchanged.

The state resulting form applying the reaction and reorientation operator <math>\mathbf{s}^{\mathcal{O}\circ\mathcal{A}}</math> is known as the post-interaction configuration and denoted by <math>\mathbf{s}^{\mathcal{I}}:=\mathbf{s}^{\mathcal{O}\circ\mathcal{A}}</math>.
[[File:Hexdynamics.png|center|thumb|888x888px|Dynamics of the BIO-LGCA model. Every time step, the occupation numbers are changed stochastically by the reaction and/or reorientation operators in all lattice sites simultaneously during the interaction step. Subsequently, particles are deterministically moved to the same velocity channel on a neighboring node in the direction of their velocity channel, during the transport step. Colors in the sketch are used to track the dynamics of the particles of individual nodes. This sketch assumes a particle-conserving rule (no reaction operator).]]

==== Transport step ====
After the interaction step, the deterministic transport step is applied synchronously to all lattice sites. The transport step simulates the movement of agents according to their velocity, due to the [[Active matter|self-propulsion]] of living organisms.

During this step, the occupation numbers of post-interaction states will be defined as the new occupation states of the same channel of the neighboring lattice site in the direction of the velocity channel, i.e. <math>s_i(\mathbf{r}+\mathbf{c}_i)=s_i^{\mathcal{I}}(\mathbf{r})</math>.

A new time step begins when both interaction and transport steps have occurred. Therefore, the dynamics of the BIO-LGCA can be summarized as the stochastic [[Recurrence relation|finite-difference]] microdynamical equation <math display="block">s_i(\mathbf{r}+\mathbf{c}_i,k+1)=s_i^{\mathcal{I}}(\mathbf{r},k)</math>

== Example interaction dynamics ==
[[File:PalignLGCA.webm|thumb|A hexagonal BIO-LGCA model of polar swarming. In this model, cells preferentially change their velocities to be parallel to the neighborhood's momentum. Lattice sites are colored according to their orientation, following the [[color wheel]]. Empty sites are white. Periodic boundary conditions were used.]]
The transition probability for the reaction and/or reorientation operator must be defined to appropriately simulate the modeled system. Some elementary interactions and the corresponding transition probabilities are listed below.

=== Random walk ===
In the absence of any external or internal stimuli, cells may move randomly without any directional preference. In this case, the reorientation operator may be defined through a transition probability<math display="block">P\left(\left.\mathbf{s}\rightarrow\mathbf{s}^{\mathcal{O}}\right|\mathbf{s}_{\mathcal{N}}\right)
=\frac{\delta_{n(\mathbf{s}),n\left(\mathbf{s}^{\mathcal{O}}\right)}}{Z}</math>
[[File:BIOLGCAbarkley.webm|thumb|382x382px|A hexagonal BIO-LGCA model of excitable media. In this model, the reaction operator favors the rapid reproduction of particles within velocity channels, and the slow death of particles within rest channels. Particles in rest channels inhibit the reproduction of particles in velocity channels. The reorientation operator is the random walk operator in the text. Lattice sites are brightly colored the more motile particles are present. Resting particles are not shown. Periodic boundary conditions were used.]]
where <math>Z=\sum_{\mathbf{s}^{\mathcal{O}}}\delta_{n\left(\mathbf{s}\right),n\left(\mathbf{s}^{\mathcal{O}}\right)}</math>. Such transition probability allows any post-reorientation configuration <math>\mathbf{s}^\mathcal{O}</math> with the same number of particles as the pre-reorientation configuration <math>\mathbf{s}</math>, to be picked [[Discrete uniform distribution|uniformly]].

=== Simple birth and death process ===
If organisms reproduce and die independently of other individuals (with the exception of the finite carrying capacity), then a simple birth/death process can be simulated<ref name=":0" /> with a transition probability given by<math display="block">P\left(\left.\mathbf{s}\rightarrow\mathbf{s}^\mathcal{A}\right|\mathbf{s}_{\mathcal{N}}\right)=
\left[r_b\delta_{n\left(\mathbf{s}^{\mathcal{A}}\right),n\left(\mathbf{s}\right)+1}
+r_d\delta_{n\left(\mathbf{s}^{\mathcal{A}}\right),n\left(\mathbf{s}\right) - 1}\right]
\Theta\left[n\left(\mathbf{s}^{\mathcal{A}}\right)\right]
\Theta\left[n\left( K-\mathbf{s}^{\mathcal{A}}\right)\right]</math> where <math>r_b,r_d\in[0,1]</math>, <math>r_b+r_d\leq 1</math> are constant birth and death probabilities, respectively, <math>\delta_{i,j}</math> is the Kronecker delta which ensures only one birth/death event happens every time step, and <math>\Theta(x)</math> is the [[Heaviside step function|Heaviside function]], which makes sure particle numbers are positive and bounded by the carrying capacity <math>K</math>.
[[File:AdhLGCA.webm|thumb|A square BIO-LGCA model of cells interacting adhesively. Cells move preferentially in the direction of the cell density gradient. Lattice sites are colored with increasingly darker blue colors with increasing cell density. Empty nodes are colored white.Periodic boundary conditions are used.]]

=== Adhesive interactions ===
Cells may adhere to one another by [[cadherin]] molecules on the cell surface. Cadherin interactions allow cells to form aggregates. The formation of cell aggregates via adhesive biomolecules can be modeled<ref>{{Cite journal|last=Bussemaker|first=Harmen J.|date=1996-02-01|title=Analysis of a pattern-forming lattice-gas automaton: Mean-field theory and beyond|url=http://dx.doi.org/10.1103/physreve.53.1644|journal=Physical Review E|volume=53|issue=2|pages=1644–1661|doi=10.1103/physreve.53.1644|pmid=9964425|bibcode=1996PhRvE..53.1644B |issn=1063-651X}}</ref> by a reorientation operator with transition probabilities defined as<math display="block">P\left(\left.\mathbf{s}\rightarrow\mathbf{s}^{\mathcal{O}}\right|\mathbf{s}_{\mathcal{N}}\right)
=\frac{1}{Z}\exp\left[\beta\mathbf{G}\left(\mathbf{s}_{\mathcal{N}}\right)\cdot\mathbf{J}\left(\mathbf{s}^{\mathcal{O}}\right)\right]</math>
[[File:ChemoLGCA.webm|thumb|A square BIO-LGCA model of cells indirectly interacting [[Chemotaxis|chemotactically]]. In this model, cells produce a diffusing chemoattractant with a certain [[Biological half-life|half-life]]. Cells preferentially move in the direction of the chemoattractant gradient. Lattice sites are additively colored with a darker blue tint with increasing cell density, and with a darker yellow tint with increasing chemoattractant concentration. Empty lattice sites are colored white. Periodic boundary conditions were used.]]
where <math>\mathbf{G}\left(\mathbf{s}_{\mathcal{N}}\right)</math> is a vector pointing in the direction of maximum cell density,  defined as <math>\mathbf{G}\left(\mathbf{s}_{\mathcal{N}}\right)=
\sum_{\mathbf{r}'\in\mathcal{N}}\left(\mathbf{r}'-\mathbf{r}\right)n\left(\mathbf{s}_{\mathcal{N}}^{\mathbf{r}'}\right)</math>, where <math>\mathbf{s}_{\mathcal{N}}^{\mathbf{r}'}</math> is the configuration of the lattice site <math>\mathbf{r}'</math> within the neighborhood <math>\mathcal{N}</math>, and <math>\mathbf{J}\left(\mathbf{s}^{\mathcal{O}}\right)</math> is the momentum of the post-reorientation configuration, defined as <math>\mathbf{J}\left(\mathbf{s}^{\mathcal{O}}\right)=\sum_{j=1}^bs_j^{\mathcal{O}}\mathbf{c}_j</math>. This transition probability favors post-reorientation configurations with cells moving towards the cell density gradient.

== Mathematical analysis ==
Since an exact treatment of a stochastic agent-based model quickly becomes unfeasible due to high-order [[Correlation and dependence|correlations]] between all agents,<ref>{{Cite journal|last1=Ovaskainen|first1=Otso|last2=Somervuo|first2=Panu|last3=Finkelshtein|first3=Dmitri|date=2020-10-28|title=A general mathematical method for predicting spatio-temporal correlations emerging from agent-based models|journal=Journal of the Royal Society Interface|volume=17|issue=171|pages=20200655|doi=10.1098/rsif.2020.0655|pmc=7653394|pmid=33109018}}</ref> the general method of analyzing a BIO-LGCA model is to cast it into an approximate, deterministic [[Recurrence relation|finite difference equation]] (FDE) describing the [[Expected value|mean]] dynamics of the population, then performing the mathematical analysis of this approximate model, and comparing the results to the original BIO-LGCA model.

First, the expected value of the microdynamical equation <math>s_m(\mathbf{r}+\mathbf{c}_m,k+1)=s_m^{\mathcal{I}}(\mathbf{r},k)</math> is obtained<math display="block">f_m\left(\mathbf{r}+\mathbf{c}_m,k+1\right)=
\left\langle s_m^{\mathcal{I}}\left(\mathbf{r},k\right)\right\rangle</math> where <math>\langle\cdot\rangle</math> denotes the [[expected value]], and <math>f_m\left(\mathbf{r},k\right):=\left\langle s_m\left(\mathbf{r},k\right)\right\rangle</math> is the expected value of the <math>m</math>-th channel occupation number of the lattice site at <math>\mathbf{r}</math> at time step <math>k</math>. However, the term on the right, <math>\left\langle s_m^{\mathcal{I}}\left(\mathbf{r},k\right)\right\rangle</math> is highly nonlinear on the occupation numbers of both the lattice site <math>\mathbf{r}</math> and the lattice sites within the interaction neighborhood <math>\mathcal{N}</math>, due to the form of the transition probability <math>P\left(\left.\mathbf{s}\rightarrow\mathbf{s}^{\mathcal{I}}\right|\mathbf{s}_{\mathcal{N}}\right)</math> and the statistics of particle placement within velocity channels (for example, arising from an exclusion principle imposed on channel occupations). This non-linearity would result in high-order correlations and moments among all channel occupations involved. Instead, a mean-field approximation is usually assumed, wherein all correlations and high order moments are neglected, such that direct particle-particle interactions are substituted by interactions with the respective expected values. In other words, if <math>X_1,X_2,\ldots,X_n</math> are random variables, and <math>F:\mathbb{R}^n\mapsto\mathbb{R}</math> is a function, then<math>\left\langle F\left(X_1,X_2,\ldots,X_n\right)\right\rangle\approx
F\left(\left\langle X_1\right\rangle,\left\langle X_2\right\rangle,\ldots,\left\langle X_n\right\rangle\right)</math> under this approximation. Thus, we can simplify the equation to<math display="block">f_m\left(\mathbf{r}+\mathbf{c}_m,k+1\right)=
\mathcal{C}\left(\mathbf{f}\left(\mathbf{r},k\right),\mathbf{f}_{\mathcal{N}}\left(\mathbf{r},k\right)\right)</math> where <math>\mathcal{C}\left(\mathbf{f}\left(\mathbf{r},k\right),\mathbf{f}_{\mathcal{N}}\left(\mathbf{r},k\right)\right)</math> is a nonlinear function of the expected lattice site configuration <math>\mathbf{f}\left(\mathbf{r},k\right)</math> and the expected neighborhood configuration <math>\mathbf{f}_{\mathcal{N}}\left(\mathbf{r},k\right)</math> dependent on the transition probabilities and in-node particle statistics.

From this nonlinear FDE, one may identify several homogeneous [[Equilibrium point|steady states]], or constants <math>\bar{f}_m</math> independent of <math>\mathbf{r}</math> and <math>k</math> which are solutions to the FDE. To study the stability conditions of these steady states and the pattern formation potential of the model, a [[Linear stability|linear stability analysis]] can be performed. To do so, the nonlinear FDE is linearized as<math display="block">f_m\left(\mathbf{r}+\mathbf{c}_m,k+1\right)=
\sum_{j=1}^K\left.\frac{\partial\mathcal{C}}{\partial f_j\left(\mathbf{r},k\right)}\right|_{\mathrm{ss}}f_j\left(\mathbf{r},k\right)+
\sum_{j=1}^K\sum_{p=1}^K\left.\frac{\partial\mathcal{C}}{\partial f_j\left(\mathbf{r}+\mathbf{c}_p,k\right)}\right|_{\mathrm{ss}}f_j\left(\mathbf{r}+\mathbf{c}_p,k\right)</math> where <math>\mathrm{ss}</math> denotes the homogeneous steady state <math>f_m\left(\mathbf{r},k\right)=\bar{f}_m,m\in\{1,\ldots,K\}</math>, and a [[von Neumann neighborhood]] was assumed. In order to cast it into a more familiar finite difference equation with temporal increments only, a [[discrete Fourier transform]] can be applied on both sides of the equation. After applying the [[Discrete Fourier transform#Shift theorem|shift theorem]] and isolating the term with a temporal increment on the left, one obtains the lattice-Boltzmann equation<ref name=":1" /><math display="block">\hat{f}_m\left(\mathbf{q},k+1\right)=e^{-\frac{2 \pi i}{L}\mathbf{q}\cdot\mathbf{c}_m}
\left\{\sum_{j=1}^K\left[\left.\frac{\partial\mathcal{C}}{\partial f_j\left(\mathbf{r},k\right)}\right|_{\mathrm{ss}}+\sum_{p=1}^K\left.\frac{\partial\mathcal{C}}{\partial f_j\left(\mathbf{r}+\mathbf{c}_p,k\right)}\right|_{\mathrm{ss}}e^{\frac{2\pi i}{L}\mathbf{q}\cdot\mathbf{c}_p}\right]\hat{f}_j\left(\mathbf{q},k\right)\right\}</math> where <math>i=\sqrt{-1}</math> is the [[imaginary unit]], <math>L</math> is the size of the lattice along one dimension, <math>\mathbf{q}\in\{1,2,\ldots,L\}^d</math> is the Fourier [[Wavenumber|wave number]], and <math>\hat{\cdot}=\mathcal{F}\{\cdot\}</math> denotes the discrete Fourier transform. In matrix notation, this equation is simplified to <math>\hat{\mathbf{f}}\left(\mathbf{q},k+1\right)=\Gamma\hat{\mathbf{f}}\left(\mathbf{q},k\right)</math>, where the matrix <math>\Gamma</math> is called the ''Boltzmann propagator'' and is defined as<math display="block">\Gamma_{m,j}=e^{-\frac{2\pi i}{L}\mathbf{q}\cdot\mathbf{c}_m}
\left[\left.\frac{\partial\mathcal{C}}{\partial f_j\left(\mathbf{r},k\right)}\right|_{\mathrm{ss}}+\sum_{p=1}^K\left.\frac{\partial\mathcal{C}}{\partial f_j\left(\mathbf{r}+\mathbf{c}_p,k\right)}\right|_{\mathrm{ss}}e^{\frac{2\pi i}{L}\mathbf{q}\cdot \mathbf{c}_p}\right].</math> The [[Eigenvalues and eigenvectors|eigenvalues]] <math>\lambda\left(\mathbf{q}\right)</math> of the Boltzmann propagator dictate the stability properties of the steady state:<ref name=":1" />

* If <math>\left|\lambda\left(\mathbf{q}\right)\right|>1</math>, where <math>|\cdot|</math> denotes the [[Complex number#Polar complex plane|modulus]], then perturbations with wave number <math>\mathbf{q}</math> grow with time. If <math>\left|\lambda\left(\mathbf{q}_{\mathrm{max}}\right)\right|>1</math>, and <math>\left|\lambda\left(\mathbf{q}_{\mathrm{max}}\right)\right|\geq\left|\lambda\left(\mathbf{q}\right)\right|\forall\mathbf{q}\in\{1,2,\ldots,L\}^d</math>, then perturbations with wave number <math>\mathbf{q}_{\mathrm{max}}</math> will dominate and patterns with a clear [[wavelength]] will be observed. Otherwise, the steady state is stable and any perturbations will decay.
* If <math>\mathrm{arg}\left[\lambda\left(q\right)\right]\neq 0</math>, where <math>\mathrm{arg}(\cdot)</math> denotes the [[Complex number#Polar complex plane|argument]], then perturbations are transported and non-stationary population behaviors are observed. Otherwise, the population will appear static at the macroscopic level.

== Applications ==
Constructing a BIO-LGCA for the study of biological phenomena mainly involves defining appropriate transition probabilities for the interaction operator, though precise definitions of the state space (to consider several cellular [[phenotype]]s, for example), boundary conditions (for modeling phenomena in confined conditions), neighborhood (to match experimental interaction ranges quantitatively), and carrying capacity (to simulate crowding effects for given cell sizes) may be important for specific applications. While the distribution of the reorientation operator can be obtained through the aforementioned statistical and biophysical methods, the distribution of the reaction operators can be estimated from the statistics of ''in vitro'' experiments, for example.<ref>{{Cite journal|last1=Dirkse|first1=Anne|last2=Golebiewska|first2=Anna|last3=Buder|first3=Thomas|last4=Nazarov|first4=Petr V.|last5=Muller|first5=Arnaud|last6=Poovathingal|first6=Suresh|last7=Brons|first7=Nicolaas H. C.|last8=Leite|first8=Sonia|last9=Sauvageot|first9=Nicolas|last10=Sarkisjan|first10=Dzjemma|last11=Seyfrid|first11=Mathieu|date=2019-04-16|title=Stem cell-associated heterogeneity in Glioblastoma results from intrinsic tumor plasticity shaped by the microenvironment|journal=Nature Communications|language=en|volume=10|issue=1|pages=1787|doi=10.1038/s41467-019-09853-z|issn=2041-1723|pmc=6467886|pmid=30992437|bibcode=2019NatCo..10.1787D }}</ref>

BIO-LGCA models have been used to study several cellular, biophysical and medical phenomena. Some examples include:

* [[Angiogenesis]]:<ref>{{Cite journal|last1=Mente|first1=Carsten|last2=Prade|first2=Ina|last3=Brusch|first3=Lutz|last4=Breier|first4=Georg|last5=Deutsch|first5=Andreas|date=2010-10-01|title=Parameter estimation with a novel gradient-based optimization method for biological lattice-gas cellular automaton models|url=http://link.springer.com/10.1007/s00285-010-0366-4|journal=Journal of Mathematical Biology|language=en|volume=63|issue=1|pages=173–200|doi=10.1007/s00285-010-0366-4|pmid=20886214|s2cid=12404555|issn=0303-6812}}</ref> an in vitro experiment with endothelial cells and BIO-LGCA simulation observables were compared to determine the processes involved during angiogenesis and their weight. They found that adhesion, alignment, contact guidance, and [[Extracellular matrix|ECM]] remodeling are all involved in angiogenesis, while long-range interactions are not vital to the process.
* [[Active fluid]]s:<ref>{{Cite journal|last1=Bussemaker|first1=Harmen J.|last2=Deutsch|first2=Andreas|last3=Geigant|first3=Edith|date=1997-06-30|title=Mean-Field Analysis of a Dynamical Phase Transition in a Cellular Automaton Model for Collective Motion|url=https://link.aps.org/doi/10.1103/PhysRevLett.78.5018|journal=Physical Review Letters|language=en|volume=78|issue=26|pages=5018–5021|doi=10.1103/PhysRevLett.78.5018|arxiv=physics/9706008|bibcode=1997PhRvL..78.5018B |s2cid=45979152|issn=0031-9007}}</ref> the macroscopic physical properties of a population of particles interacting through polar alignment interactions were investigated using a BIO-LGCA model. It was found that increasing initial particle density and interaction strength result in a [[Second-order phase transition|second order phase transition]] from a homogeneous, disordered state to an ordered, patterned, moving state.
* [[Epidemiology]]:<ref name=":2">{{Cite journal|last1=Fuks|first1=Henryk|last2=Lawniczak|first2=Anna T.|date=2001|title=Individual-based lattice model for spatial spread of epidemics|journal=Discrete Dynamics in Nature and Society|volume=6|issue=3|pages=191–200|language=en|doi=10.1155/s1026022601000206|doi-access=free|hdl=1807/82157|hdl-access=free}}</ref> a spatial [[SIR Model|SIR]] BIO-LGCA model was used to study the effect of different vaccination strategies, and the effect of approximating a spatial epidemic with a non-spatial model. They found that barrier-type vaccination strategies are much more effective than spatially uniform vaccination strategies. Furthermore, they found that non-spatial models greatly overestimate the rate of infection.
* Cell [[Jamming (physics)|jamming]]:<ref>{{Cite journal|last1=Ilina|first1=Olga|last2=Gritsenko|first2=Pavlo G.|last3=Syga|first3=Simon|last4=Lippoldt|first4=Jürgen|last5=La Porta|first5=Caterina A. M.|last6=Chepizhko|first6=Oleksandr|last7=Grosser|first7=Steffen|last8=Vullings|first8=Manon|last9=Bakker|first9=Gert-Jan|last10=Starruß|first10=Jörn|last11=Bult|first11=Peter|date=2020-08-24|title=Cell–cell adhesion and 3D matrix confinement determine jamming transitions in breast cancer invasion|journal=Nature Cell Biology|language=en|volume=22|issue=9|pages=1103–1115|doi=10.1038/s41556-020-0552-6|issn=1476-4679|pmc=7502685|pmid=32839548}}</ref> in vitro and Bio-LGCA models were used for studying [[Metastasis|metastatic]] behavior in breast cancer. The BIO-LGCA model revealed that metastasis may exhibit different behaviors, such as random gas-like, jammed solid-like, and correlated fluid-like states, depending on the adhesivity level among cells, ECM density, and cell-ECM interactions.

== References ==

{{reflist}}

== External links ==
* [https://imc.zih.tu-dresden.de//biolgca/ Bio-LGCA Simulator] - An online simulator with elementary interactions with personalizable parameter values.
* [https://github.com/sisyga/BIO-LGCA BIO-LGCA Python Package] - An open source Python package for implementing BIO-LGCA model simulations.

[[Category:Statistical mechanics]]
[[Category:Lattice models]]
[[Category:Stochastic models]]
[[Category:Complex dynamics]]