*Editorial* **Applying the Free Energy Principle to Complex Adaptive Systems**

**Paul B. Badcock 1,2,\*, Maxwell J. D. Ramstead 3,4, Zahra Sheikhbahaee <sup>5</sup> and Axel Constant <sup>6</sup>**


The free energy principle (FEP) is a formulation of the adaptive, belief-driven behaviour of self-organizing systems that gained prominence in the early 2000s as a unified model of the brain [1,2]. Since then, the theory has been applied to a wide range of biotic phenomena, extending from single cells and flora [3,4], the emergence of life and evolutionary dynamics [5,6], and to the biosphere itself [7]. For our part, we have previously proposed that the FEP can be integrated with Tinbergen's seminal four questions in biology to furnish a multiscale ontology of living systems [8]. We have also explored more specific applications, e.g., to the evolution and development of human phenotypes [9–11], sociocultural cognition, behaviour, and learning [12,13], as well as the dynamic construction of environmental niches by their denizens [14,15].

Despite such contributions, the capacity of the FEP to extend beyond the human brain and behaviour, and to explain living systems more generally, has only begun to be explored. This begs the following questions: Can the FEP be applied to any organism? Does it allow us to explain the dynamics of all living systems, including large-scale social behaviour? Does the FEP provide a formal, empirically tractable theory of any complex adaptive system, living or not? With such questions in mind, the aim of this Special Issue was to showcase the breadth of the FEP as a unified theory of complex adaptive systems, biological or otherwise. Instead of concentrating on the human brain and behaviour, we welcomed contributions that applied the FEP to other complex adaptive systems, with the hope of exemplifying the extent of its explanatory scope.

For the uninitiated, it is worth briefly outlining what the FEP is. Variational free energy refers to an information theoretic quantity that places an upper limit on the entropy of a system's observations, relative to a generative model instantiated by an agent. (In this context, entropy is defined as the time-average of 'surprise' or the negative log probability of the agent's sensory data.) Generative models harness probabilistic mappings from hidden causes in the environment to observed consequences (i.e., sensory data), and state transitions inherent to the environment [2]. Under the FEP, an organism is modelled as implicitly 'expecting' to find itself within a limited range of phenotypic states; as such, deviations from these states elicit a type of 'phenotypic surprise' (i.e., the deviation between actual and phenotypically preferred states). Consequently, organisms remain alive by acting in ways that minimize this type of surprise (e.g., a fish avoiding the 'state' of being out of water). In other words, and more heuristically, free energy scores the discrepancy between desired and sensed data; and the FEP states that the imperative of all self-organizing

**Citation:** Badcock, P.B.; Ramstead, M.J.D.; Sheikhbahaee, Z.; Constant, A. Applying the Free Energy Principle to Complex Adaptive Systems. *Entropy* **2022**, *24*, 689. https:// doi.org/10.3390/e24050689

Received: 6 May 2022 Accepted: 11 May 2022 Published: 13 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

<sup>1</sup> Centre for Youth Mental Health, The University of Melbourne, Melbourne, VIC 3010, Australia

systems is to keep this discrepancy at bay by bringing about preferred observations via action (see Figure 1).

**Figure 1. The free energy principle**. (**A**) Schematic of the quantities that characterise free energy, including a system's internal states, *μ* (e.g., a brain), and quantities that describe the system's exchanges with the environment; specifically, its sensory input, *s* = *g*(*η*,*a*) + *ω*, and actions, *a*, which alter the ways in which the system samples its environment. Environmental states are further defined by equations of motion, . *η = f(η,a) + ω*, which describe the dynamics of (hidden) states extraneous to the system, *η*, whereas *ω* refers to random fluctuations. Under this scheme, internal states and action operate synergistically to reduce free energy, which reflects a function of sensory input and the probabilistic representation (variational density), *q*(*η*:*μ*), that internal states encode. Note that external and internal states are statistically separated by a Markov blanket, which possesses both 'sensory' and 'active' states. Internal states are influenced by, but cannot affect, sensory states, whereas external states are influenced by, but cannot affect, active states, creating a conditional independence between the system and its environment. (**B**) Alternative equations that describe the minimisation of free energy. With respect to action, free energy can only be suppressed by the system's selective sampling of (predicted) sensory input, which increases the accuracy of its predictions. On the other hand, optimising internal states minimises divergence by making the representation an approximate conditional density on the hidden causes of sensory input. This optimisation reduces the free energy bound on surprise, which means that action allows the system to avoid surprising sensations. Reproduced from [8].

There are two main ways for a self-organizing system to minimize free energy. The first is by changing its perception of the world. Previously, this has been explored through reference to human neural processing. The FEP appeals to a view of the brain as a hierarchical "inference machine", which instantiates a hierarchy of hypotheses about the world (i.e., a generative model) that enables an organism to minimize free energy (and therefore keep entropic dissipation at bay, at least locally) by reducing discrepancies between incoming sensory inputs and top-down predictions (i.e., prediction errors) [2]. Neurobiologically xpectations about sensory data are thought to be encoded by deep pyramidal cells (i.e., representation units) at every level of the cortical hierarchy, which carry predictions downward to suppress errors at the level below, whereas prediction errors themselves are encoded by superficial pyramidal cells and are carried forward to revise expectations at the level above [16]. The relative influence of ascending (error) and descending (representation) signals is weighted by their inverse variance or *precision* (e.g., a high precision on ascending error signals lowers confidence in descending predictions), which is mediated by neuromodulation and reflected psychologically by attentional selection and sensory attenuation. In short, the recursive neural dynamics described here enable us to minimise free energy (resp. prediction error) by updating our internal models (i.e., perception).

Second, an organism can reduce surprise directly by acting upon the world in order to fulfill its expectations and generate unsurprising sensations. This process of 'active inference' describes how an organism reduces free energy through self-fulfilling cycles of action and perception [17]. Active inference models implement action selection as the minimization of *expected free energy*, which is the free energy expected under beliefs about possible courses of action, or policies. By selecting actions that are expected to minimize free energy, the organism can maintain itself within preferred, phenotypically unsurprising states. Thus, survival mandates that an organism must not only reduce free energy from moment to moment; it must also reduce the expected free energy associated with the future outcomes of action [18,19].

Having briefly outlined the rudiments of the FEP, let us turn briefly to complex adaptive systems (CAS). This concept is synonymous with complexity science and has its roots in evolutionary systems theory, which assumes a dynamic, inextricable relationship between generalised selection and self-organization [11]. Broadly speaking, a CAS refers to any multi-component, self-organising system that adapts to it environment through an autonomous process of selection, which recruits the outcomes of localised interactions between its components to select a subset of those components for replication and enhancement [20]. Holland [21] describes four key features of CAS: they consist of large numbers of components that interact by sending and receiving signals (i.e., *parallelism*); the actions of their components depend upon the signals they receive (i.e., *conditional action*); groups of rules can form subroutines that can be combined to deal with environmental novelties (i.e., *modularity*); and the components of the system change over time to adapt to the environment and improve performance (i.e., *adaptation*). Applications of the CAS framework have proliferated across the physical, human and computer sciences, but there is not the scope to survey this literature here. However, to pre-empt the papers to follow, we would note that this framework has already been applied to precisely the same systems that are the foci of our contributors–ranging from metabolic and cellular processes, e.g., [22–26]; to the brain and social processes, e.g., [26–30]; and to artificial intelligence and robotics, e.g., [31–33]. The articles presented in our Special Issue build upon such literature by illustrating how the FEP can afford fresh insights into the dynamics of CAS.

Three of the contributions to our Special Issue leverage the FEP to cast new light on processes intrinsic to biological agents. In *Cancer Niches and their Kikuchi Free Energy*, Sajid, Convertino et al. [34] examine cancer morphogenetic fields as undesirable stable attractors in the complex dynamics of homeostasis, self-renewal and differentiation, which contributes to their deviation from regular autopoietic homeostasis (the internal molecular dynamics that regulate the production and regeneration of a system's components). Sajid, Convertino et al. offer a computational model in silico to study communication and information processing at a population level of cancer cell networks within their environment in vivo. By deploying the Kikuchi free energy approximation, which is a generalisation of the Bethe free energy for computing beliefs over large sets of cell clusters, they

account for higher-order interactions and phase transitions between clusters of healthy and oncogenic cells. Here, cancer niche construction can be construed as a Bayes-optimal process for the transmission of information across different levels of cellular networks due to its tendency to minimize the overall Kikuchi free energy over the whole system. Their findings suggest that three distinct cancer trajectories–namely, proliferation (local growth), metastasis and apoptosis–can emerge from the natural evolution of the state function (i.e., free energy) in biological systems. These findings have important implications for our understanding and study of cancer cell growth and apoptosis.

Next, Parr describes how biochemical networks in adaptive biological systems can be recast in terms of an inferential message passing scheme that involves the gradient descent on variational free energy towards the least surprising states, based on the organism's implicit (generative) model of these states. In *Message Passing and Metabolism* [35], he points out that the biochemical regulation of metabolic processes relies on sparse interactions (message-passing) between coupled reactions, with enzymes creating conditional dependencies between reactants. He then extrapolates the *law of mass action* (the rate of chemical reaction and the concentrations of reactants involved in this process) and the *Michaelis–Menten kinetics* (which approximates the dynamics of irreversible enzymatic reactions) from the FEP. Assuming the existence of a causal structure in biochemical (metabolic) networks, one can build the sparse message passing scheme to capture the independence of substrates and products, conditioned upon the enzyme and enzyme-substrate complex within such networks. The temporal evolution of the categorical probability of each state within this system can be described by a chemical master equation that takes into account sparse network interactions. Parr describes how the steady state distribution of these dynamics can be recast as a generative model, which suggests that the biochemistry that underlies metabolism follows an inferential message-passing scheme that seeks to minimise free energy. An important extension of Parr's model is that metabolic disorder can emerge when an enzymatic disconnection by thiamine depletion interrupts message passing and incites aberrant prior beliefs, which gives rise to false (biochemical) inference.

The third contribution follows more traditional applications of the FEP by accounting for conscious, first-person experience. In *The Radically Embodied Conscious Cybernetic Bayesian Brain*, Safron [36] proposes models of embodied conscious agency based on the FEP, extending the Integrated World Modelling Theory of consciousness proposed elsewhere [37] to explicitly account for aspects of intentional actions and agentic experiences. According to the radically embodied account on offer, what we call attention and imagination emerge from the (sometimes liminal) activity of multimodal, action-oriented body maps and representations, realized as neural attractors in the form of 'embodied selfmodels' (ESMs), which conform to the FEP as cybernetic controllers. When functioning online, ESMs allow for overt interactions with affordances, or structured possibilities for environmental interactions. However, Safron suggests subthreshold activations of such 'quasi-homuncular' ESMs also underwrite our (affordance-structured) covert abilities to imagine and pursue courses of action, as well as our ability to intentionally deploy attentional resources. Thus, even seemingly abstract representational capacities may be grounded in twin capacities for embodied action and counterfactual explorations of the world. Safron then applies this radically embodied perspective to core aspects of conscious experience. He attempts to chart a middle way between perspectives in the representation wars in cognitive science, describing brains as hybrid machine learning architectures capable of supporting both symbolic and sub-symbolic processes for 4E agents (where cognition is thought to be *embodied*, *embedded*, *enacted* and *extended*). Safron's perspective is ecumenical, deploying information-theoretic constructs and representationalist concepts that would be rejected by hard-line proponents of both 4E cognition and more Cartesian (representationalist) approaches. For example, to account for information flow in mammalian brains, Safron deploys constructs that are typically rejected by 4E theorists, such as Cartesian theatres and quasi-homunculi. However, he does so from a radically embodied

perspective, suggesting that such a "strange inversion of reasoning" follows from principles of cognitive development and computational neuroscience.

Unlike the authors above, Goekoop and de Kleijn look beyond the phenotype to consider how the FEP might apply to groups. In *Permutation Entropy as a Universal Disorder Criterion* [38], they argue that living systems can be described as hierarchical problemsolving machines that embody predictive models of their econiches, called a goal hierarchy, which incorporates a set of lower-order econiches (goals) and corresponding subniches (subgoals) that the system needs to pursue in order to achieve the global econiche (goal) represented at the top of the hierarchy. Using this scheme to frame the rest of their argument, they concentrate on stress responses in organisms, dyads and collectives. Equating stress with free energy or 'prediction error', and stress responses with 'action', they suggest that as free energy increases, there is a progressive collapse of (allostatic) hierarchical control, eventually resulting in disordered states characterised by behavioural shifts from longterm goal-directed behaviour (e.g., reproductive success) towards short-term goals and habitual behaviours concentrated on self-preservation (e.g., survival). After introducing permutation entropy as a universal measure of disordered states across such systems, they briefly describe how their model can be used to explain disorder at an individual level, before progressing to the transmission of disorder through interpersonal interactions, and concluding with a brief discussion of population-level dynamics.

The idea that the FEP can be extended to social systems is also taken up by Kaufmann and colleagues. In *An Active Inference Model of Collective Intelligence* [39], the authors propose an active inference model of alignment, describing the manner in which withinscale local interactions (e.g., individual agents' behaviors) can align with cross-scale global phenomena (e.g., collective behavior) in multi-scale systems. In so doing, they offer a principled, agent-based model that has the potential to function as a workbench to simulate collective intelligence as an emergent phenomenon, across many scales. Although one obvious target for this modelling approach would be human behavior as an emergent phenomenon that ties physiological, cognitive and cultural dynamics, nothing, in principle, limits the application of Kaufmann and colleagues' model to human phenomena.

In another paper that illustrates the broad scope of multiscale thinking under active inference, Jesse Hoey, in *Equality and Freedom as Uncertainty in Groups* [40], shows how agents attempting to align with other group members leads to a quasi-equilibrium, or "sweet spot", at which the group free energy is minimal and the agent's predictive capacity of higher order parameters, such as those attributed to the group, matches the group's capacity to predict an agent's behavior. Hoey further discusses two intriguing trade-offs. Higher agent model complexity leads to lower individual learning capacity with respect to the complexity of the group, resulting in agents who are hobbled in the pursuit of their own ends, but in a group that is more diverse, innovative, and open to change. On the other hand, lower agent model complexity allows the expression of individual preference towards the group, but the group becomes more homogeneous, secure, and closed, as otherwise the pro-social behaviours of individual agents would be hampered. Hoey suggests that such emergent social dynamics provides insights into concepts such as freedom and equality in society, which correspond to changes in model uncertainties and complexities. Oscillation between radical freedom, where no cooperation is possible, and radical equality, where no discovery is possible, is an emergent phenomenon characteristic of Western society; akin to what Karl Marx called historical materialism, which according to many is the main driver of history itself. Could Hoey- s findings initiate research on an active inference model of history as an emergent phenomenon of human societies?

So far, we have considered a range of applications of the FEP to living systems. However, active inference–the process theory derived from the FEP–is increasingly being applied to machine intelligence in practical settings. Use cases in robotics provide an exciting opportunity to test the applicability of active inference to implement sensory processes and motor control in real time. In their review for our Special Issue, *How Active Inference Could Help Revolutionise Robotics*, Da Costa and colleagues [41] examine the

usefulness of active inference for several core problems in robotics, such as state estimation in artificial perception, motor control, learning, safety and explainability. They argue that active inference may help advance robotics due to several of its core features: it enables explainable artificial intelligence in a manner that operates in situations involving uncertainty, volatility, and noise. This is especially relevant for human-centric applications, such as human–robot interaction and healthcare.

In closing, it is worth recognising that the majority of submissions presented herein focus chiefly on human systems, despite our call for more wide-ranging applications. Nevertheless, it should be clear that the authors' proposals can be readily extended to other complex adaptive systems, including biological dynamics intrinsic to other lifeforms [34–36], collective, group-level behaviour [39,40], and even non-living systems [41]. Taken together, we hope that the collection of papers presented in our Special Issue motivate others to consider how the FEP might be gainfully applied to their own systems of interest, living or otherwise.

**Funding:** ZS is supported by The Long-Term Future Fund, provided by the Centre for Effective Altruism. AC is supported by the Australian Laureate Fellowship project, A Philosophy of Medicine for the 21st Century (Ref: FL170100160), and by a Social Sciences and Humanities Research Council doctoral fellowship (Ref: 752–2019-0065).

**Acknowledgments:** We would like to thank Lancelot Da Costa, Rutger Goekoop, Jesse Hoey, Rafael Kaufmann, Thomas Parr, Adam Safron, and Noor Sajid for their comments on an earlier draft of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Cancer Niches and Their Kikuchi Free Energy**

**Noor Sajid 1,\*,†, Laura Convertino 1,2,\*,† and Karl Friston <sup>1</sup>**


**Abstract:** Biological forms depend on a progressive specialization of pluripotent stem cells. The differentiation of these cells in their spatial and functional environment defines the organism itself; however, cellular mutations may disrupt the mutual balance between a cell and its niche, where cell proliferation and specialization are released from their autopoietic homeostasis. This induces the construction of cancer niches and maintains their survival. In this paper, we characterise cancer niche construction as a direct consequence of interactions between clusters of cancer and healthy cells. Explicitly, we evaluate these higher-order interactions between niches of cancer and healthy cells using Kikuchi approximations to the free energy. Kikuchi's free energy is measured in terms of changes to the sum of energies of baseline clusters of cells (or nodes) minus the energies of overcounted cluster intersections (and interactions of interactions, etc.). We posit that these changes in energy node clusters correspond to a long-term reduction in the complexity of the system conducive to cancer niche survival. We validate this formulation through numerical simulations of apoptosis, local cancer growth, and metastasis, and highlight its implications for a computational understanding of the etiopathology of cancer.

**Keywords:** cancer niches; free energy; Kikuchi approximations; apoptosis; metastasis; cluster variation method

#### **1. Introduction**

The human body develops via the progressive specialisation of pluripotent stem cells, niche construction, and organ development (i.e., morphogenesis). Understanding these processes is not only fundamental to understanding how multicellular organisms come to be, but also has important implications for aberrant events. Cellular de-differentiation can cause (and be a consequence of) disruption of the mutual balance between cells and their niche, where cell proliferation and genetic regulation are released from their autopoietic homeostasis. The disruption of this balance and the creation of new niches drives cancer cell growth and sustains their survival. The concept of a stem cell niche [1,2] identifies the "whole" of the wholeness of the microenvironment in which stem cells differentiate, grow, and survive, and the humoral, paracrine, physical, metabolic, neuronal, structural properties with which the cell exchanges information. These specialised microenvironments are found in adult organisms and their homeostatic disruption may facilitate the development of cancer colonies through an interaction between cells and their niches [3]. In this paper, we pursue the notion that cancer niche construction requires individual cells with oncogenic potential to interact in subpopulations (i.e., clusters) of healthy and cancerous cells to reach a particular attractor state.

Briefly, cancer cell niche construction results from a maladaptive degeneration of the complex dynamic homeostasis that undergirds the balanced survival, renewal, and repair of mature organs. This homeostasis is replaced by a new niche that is beneficial to cancer growth, without which oncogenic cells might not be able to escape the physiological

**Citation:** Sajid, N.; Convertino, L.; Friston, K. Cancer Niches and Their Kikuchi Free Energy. *Entropy* **2021**, *23*, 609. https://doi.org/10.3390/ e23050609

Academic Editor: Geert Verdoolaege

Received: 30 March 2021 Accepted: 7 May 2021 Published: 14 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

mechanisms of control and elimination of the degenerated cells [4]. We acknowledge the phenotypical differences between cancer cells and cancer stem cells, but for simplicity, we operationalise all cell types as simply cancer cells. Note, that cancer stem cells are a small subpopulation of cancer cells with capability of self-renewal, differentiation, and tumorigenesis [4,5].

To become carcinogenic, cells acquire some key properties that free them from their physiological life cycle. The main properties are: (i) self-sufficiency in growth signals; (ii) insensitivity to anti-growth signals; (iii) tissue invasions and metastasis; (iv) limitless reproductive potential; (v) angiogenesis; and (vi) avoiding apoptosis [6]. It may not be necessary for a cell to acquire all these properties to become a cancer cell, since different combinations of these properties are sufficient for inducing oncogenesis. For example, only some cancers have metastatic properties, namely the ability to leave the tissue of origin, move to another location and colonise another organ. The presence of metastatic sites determines the progression of cancer, staging, and patient prognosis. Fortunately, healthy organisms can fight oncogenesis via different mechanisms of self-preservation. Among these, apoptosis (i.e., programmed death) is crucial. This is a controlled mechanism of cell death that intervenes in normal organ repair, growth, renewal, as well as in inflammation and cancer.

For a computational understanding of cancer niche construction, a 2D Ising model has been used to evaluate phase transitions between cancerous and healthy cells [7–9]. For example, [10] used the Ising Hamiltonian to model metastasis, where the updates in the energy function were modelled using mean-field approximations, i.e., only considering average interactions. This factorised approach is consistent with other applications of the Ising model for simulating cancer progression [7–9]; however, cell niches are the result of numerous mechanisms, both at the cellular and population level. The single cells are involved in the construction and maintenance of the niche and they adapt via genetic, epigenetic, metabolic, structural, internal, and external signalling mechanisms. More recently, theoretical work has underlined the fundamental role of population dynamic and mutual information exchange in guiding the fate of local subpopulations of cells and the niche as a whole [11–13] via cell-to-cell, cell-to-niche, and niche-to-cell information flow.

Consequently, our work builds upon previous computational approaches by (i) casting cancer niche construction as a direct consequence of interactions between clusters of cancer and healthy cells, and (ii) using a Kikuchi Hamiltonian as a way to account for higher-order interactions when evaluating the state functions of such systems [14–19]. This allows us to move away from thinking about cells in isolation and towards accounting for interactions within cancer niches. Briefly, Kikuchi formulation approximates the free energy as a sum of local free energies of a cluster of cells (or nodes) in the system. In doing so, it provides a way to define a local population (or base cluster) which includes all the interactions between the cells.

In what follows, we simulate three cancer trajectories to provide proof of the principle of a computational formulation for oncogenic etiopathology. For this, we simulate a 2D Ising model for evaluating how cancer trajectories unfold using Kikuchi free energy approximation. Explicitly, we manipulate four (hyper)parameters, namely an interaction parameter that regulates the type of cell interaction, a tolerance parameter that determines the acceptance level of cancerous cells, a growth parameter that regulates the number of cell states switched during a single trial, and a noise parameter that influences the transition dynamics. We describe the technical details in Section 3.

First, we simulate local cancer growth within the tissue of origin by allowing pairwise interactions between cancerous and healthy cells (i.e., a high interaction parameter value). This allows us to reproduce the homeostatic shift in a healthy organ that results from the acquired oncogenic properties of single cells within it, but it is crucially sustained and enabled by the broader dynamic at the sub-population level. Secondly, we model the metastasis where cells exit their primary site and invade secondary locations. Here, cells have to prepare for the new environment before arrival in the metastatic site. Without a favourable

environment, the cells would not be able to colonise the new organ. We treat this cancer spread as an embedded mechanism in the niche construction process. This is achieved by allowing pairwise interactions between/across cancerous and healthy cells (i.e., low positive interaction parameter value), alongside a restricted capacity to grow outside the initial cancer site (i.e., decaying noise parameter). Finally, we consider apoptosis, namely the programmed death of the cancer cells. In our simulations, we assume that the cancer cells do not acquire the ability to escape apoptosis, and the physiological mechanisms of control activate to eliminate the pathological cells (via a low tolerance threshold).

The remainder of the paper is organised as follows. In Section 2, we briefly review the notion of free energies, specifically their relevance in evaluating systems under thermodynamic equilibrium. We then describe and provide simple examples for the use of Kikuchi approximations in evaluating a system state function. Equipped with this, Section 3 describes the system used to simulate cancer niche construction, and the parameters required to simulate cancer progression. Section 4 details the simulation results for the three cancer trajectories. We conclude by highlighting the implications for specific etiopathology, and potential future directions.

#### **2. Free Energies**

In this section, we review the notion of free energy and its relevance for modelling cancer niche construction. For this, we follow the formulation introduced in [20]. Briefly, free energy is the state function of a (random dynamical) system that possesses a steadystate or pullback attractor. Its value is determined by the current state of the system. The free energy can be used to describe the spatiotemporal evolution of a system as it converges to an equilibrium or nonequilibrium steady state, e.g., a particular tissue population in the human body or the body itself. The distinction between nonequilibrium and equilibrium steady state rests upon the presence of solenoidal or divergence free flow. In the absence of solenoidal dynamics, the flow of systemic states is entirely dissipative, and the steady state is at (thermodynamic) equilibrium.

To make this concrete, we introduce a closed 2D system (e.g., an Ising model) with a set of *N* discrete random variables, {*X*1, *X*2, *X*3,..., *XN*} arranged in a square grid graph (Figure 1). Each random variable, *Xi*, has a possible realisation *xi*. Here, we cast these realisations as distinct cell states in the body, or a particular tissue population, that can be either healthy or cancerous. The overall system state is denoted by the vector **x** = {*x*1, *x*2,..., *xN*}, with the corresponding energy of the system given by its Hamiltonian, *H*(**x**). For computational ease, we deal with an effective Hamiltonian to represent the system in a reduced space through nonlinear averaging of the true Hamiltonian. Consequently, this only describes a part of the eigenvalue spectrum of the true Hamiltonian. This formulation is in contrast to the molecular Hamiltonian where the Hamiltonian is decomposed into two or more separable parts (i.e., nuclei and electrons), and their interactions. For our purposes, this separation could involve the inclusion of additional random variables (that represent external states or distinct internal states) (**z**) with a Hamiltonian of the following form: *H* = *H*1(**x**) + *H*2(**z**). Practically, decompositions of the Hamiltonian of this sort are potentially important because they could naturally account for local interactions; however, their inclusion would not speak to the long-range interactions that are the current focus.

At steady state, the probability of finding the system in this state is given by Boltzmann's law, which can be expressed as [21]:

$$\begin{aligned} P(\mathbf{x}) &= \frac{1}{Z} e^{-\beta H(\mathbf{x})} \\ \ln P(\mathbf{x}) &= -\beta H(\mathbf{x}) - \ln Z \\ Z &\stackrel{\scriptstyle \Delta}{=} \sum\_{\mathbf{x} \in M} e^{-\beta H(\mathbf{x})} \end{aligned} \tag{1}$$

where *β* is the inverse temperature or precision that shapes the probability density of the distribution, *Z* is the partition function, and *M* is the set of all possible states of the system. For our purpose, we make simplifying assumptions about our system. and set *β* = 1. Briefly, we assume that the joint probability distribution describes some nonphysical system. This allows us to view Boltzmann's law as a way of defining the energy of the system where is simply an arbitrary unit that scales the energy measure [20].

**Figure 1.** Schematic illustration of the 2D system with a set of *N* discrete random variables where *N* = 1282. Here, the circles denote the variable node for each variable (*xi*), squares denote the factor nodes (*ϕa*), and edges connect the variable node *i* to the factor *a*.

The partition function, *Z*, is closely related to Helmholtz free energy, *FH* [20,22,23]:

$$\begin{split} F\_H &= -\ln Z \\ &= \ln P(\mathbf{x}) + H(\mathbf{x}) \end{split} \tag{2}$$

This quantity can be approximated using variational calculus [24], which allows an otherwise intractable *FH* to be approximated. This requires the introduction of a variational distribution, *Q*(**x**), assuming ∑ **x**∈*M Q*(**x**) = 1 and 0 ≤ *Q*(**x**) ≤ 1; ∀**x**, and the corresponding variational (or Gibbs) free energy, *FV* [20]:

$$\begin{aligned} F\_V &\equiv \underbrace{\mathcal{U}(Q)}\_{\text{Internal energy}} - \underbrace{\mathcal{S}(Q)}\_{\text{Entropy}} \\ &= \sum\_{\mathbf{x} \in \mathcal{M}} Q(\mathbf{x}) H(\mathbf{x}) - \left[ - \sum\_{\mathbf{x} \in \mathcal{M}} Q(\mathbf{x}) \ln Q(\mathbf{x}) \right] \\ &= \sum\_{\mathbf{x} \in \mathcal{M}} Q(\mathbf{x}) [\ln Q(\mathbf{x}) - \ln P(\mathbf{x})] + F\_H \\ &= \underbrace{D\_{KL} \left[ Q(\mathbf{x}) || P(\mathbf{x}) \right]}\_{\text{Complexivity}} + F\_H \end{aligned} \tag{3}$$

where *H*(**x**) = − ln *P*(**x**) + *FH*, and *DKL* is the Kullback–Leibler divergence. We know that *DKL*[*Q*(**x**)||*P*(**x**)] ≥ 0 [25], hence *FV* is an upper bound on the *FH*:

$$F\_V \ge F\_H \tag{4}$$

Thus, it follows that by minimising *FV* with respect to *Q*(**x**), we can get a good approximation of *FH* and recover *p*(**x**); however, as *N* → *c*, *c* 1, this also becomes intractable. Therefore, a practical solution is to consider the upper bound *FH* by minimising *FV* with respect to a restricted class of variational distributions *Q*(**x**). A standard restriction over the variational distribution follows from the mean-field approach, i.e., an absence of interactions [26–28]:

$$\begin{aligned} Q\_{MF}(\mathbf{x}) &= \prod\_{i=1}^{N} Q\_i(\mathbf{x}\_i) \Rightarrow \\ F\_{MF} &= \mathcal{U}\_{MF}(Q) - S\_{MF}(Q) \\ &= -\sum\_{\mathbf{z}} \sum\_{\mathbf{x}\_d} \ln \, q\_a(\mathbf{x}\_d) \prod\_{i \in N(d)} Q\_i(\mathbf{x}\_i) + \sum\_{i=1}^{N} \sum\_{\mathbf{x}\_i} Q\_i(\mathbf{x}\_i) \ln Q\_i(\mathbf{x}\_i) \\ H(\mathbf{x}) &= -\sum\_{\mathbf{z}} \ln \, q\_a(\mathbf{x}\_d) \end{aligned} \tag{5}$$

where the Hamiltonian is defined as the sum of its factors, *ϕa*, under a factor graph probability distribution function (Figure 1). Conversely, we could consider more complicated factorisations, like structure mean-field approaches [29], Bethe free energy [15,30], or Kikuchi free energy approximations or a cluster variational method [16,17], which can provide more accurate approximations. Interestingly, the Kikuchi free energy is a generalisation of the Bethe free energy, using higher-order approximations [17].

#### *Kikuchi Free Energy*

In this work, we use Kikuchi approximation to evaluate the free energy of the system [15–17,20,31]. This allows us to account for higher-order interactions between variable nodes, mimicking the type of interactions present during cell niche construction. Practically, this characterises the system evolution in terms of interactions between neighbours of cancerous and health nodes, up to size *d*. This is calculated as changes to the sum of the energies of baseline clusters of variable nodes, minus the overcounted interactions (and interactions of interactions). The premise is that in accounting for higher-order interactions within clusters, we can get a better evaluation of the systems overall state function (i.e., free energy).

Formally, the Kikuchi free energy *FK* is:

$$F\_K = \sum\_{r \in R} c\_r F\_r \tag{6}$$

where:

$$\begin{split} F\_{\mathbf{r}} &= \mathcal{U}\_{\mathbf{r}}(Q) - S\_{\mathbf{r}}(Q) \\ &= \sum\_{\mathbf{x}\_{\mathbf{r}}} Q\_{\mathbf{r}}(\mathbf{x}\_{\mathbf{r}}) H\_{\mathbf{r}}(\mathbf{x}\_{\mathbf{r}}) + \sum\_{\mathbf{x}\_{\mathbf{r}}} Q\_{\mathbf{r}}(\mathbf{x}\_{\mathbf{r}}) \ln Q\_{\mathbf{r}}(\mathbf{x}\_{\mathbf{r}}) \\ H\_{\mathbf{r}}(\mathbf{x}\_{\mathbf{r}}) &= -\sum\_{a \in A\_{\mathbf{r}}} \ln \varphi\_{a}(\mathbf{x}\_{a}). \end{split} \tag{7}$$
 
$$\mathbf{c}\_{\mathbf{r}} = 1 - \sum\_{s \in supr\mathbf{r}(\mathbf{r})} \mathbf{c}\_{\mathbf{s}} \tag{8}$$

where *r* denotes a region, i.e., the set of variable nodes within a cluster of size *d*, and *R* a finite set of all possible regions. Additionally, *cr* is the overcounting number of regions, and super(*r*) is the set of all super-regions of *r*. This follows from how *FK* is approximated. Generally, for each factor node within a cluster we need to include all adjacent variables nodes and sum the computed free energy; however, this process results in repeated counting of the interactions between sets within a cluster. Therefore, we need to subtract the free energy of interactions and interactions of interactions.

Figure 2 shows an example Kikuchi approximation using part of the system defined in Figure 1 for two different cluster sizes (2,3). To evaluate the Bethe approximation (i.e., *d* = 2) we use pairwise clusters. The Bethe entropy would then be the sum of all entropies in the pairwise cluster, subtracting the overcounted cluster interactions [15]:

$$\begin{array}{l} S\_{d=2}(\mathbf{x}) &= S(\mathbf{x}\_{1,2}) + S(\mathbf{x}\_{2,3}) + S(\mathbf{x}\_{129,130}) + S(\mathbf{x}\_{130,131}) \\ &+ S(\mathbf{x}\_{1,129}) + S(\mathbf{x}\_{2,130}) + S(\mathbf{x}\_{3,131}) \\ &- S(\mathbf{x}\_{1}) - S(\mathbf{x}\_{129}) - S(\mathbf{x}\_{3}) - S(\mathbf{x}\_{131}) \\ &- 2S(\mathbf{x}\_{2}) - 2S(\mathbf{x}\_{130}) \\ \end{array} \tag{8}$$
 
$$\mathcal{S} \quad = -\sum\_{i} Q(\mathbf{x}\_{i}) \ln Q(\mathbf{x}\_{i})$$

where *x*<sup>1</sup> appears in two clusters, {*x*1,2}, {*x*1,129}, i.e., its entropy was overcounted only once and we subtract it once. Conversely, *x*<sup>2</sup> appears in three clusters, {*x*1,2}, {*x*2,3}, {*x*2,130}, and this mean that its entropy was overcounted twice, and we subtract it twice.

*xi*

**Figure 2.** Example of Kikuchi free energy with two different cluster sizes. Panel **A** is for *d* = 2 or Bethe approximation and panel **B** is for *d* = 3.

We could approximate *FK* differently by using a cluster size of 3:

$$\begin{aligned} S\_{d=3}(\mathbf{x}) &= S(\mathbf{x}\_{129,1,2}) + S(\mathbf{x}\_{129,130,2}) + S(\mathbf{x}\_{130,2,3}) + S(\mathbf{x}\_{130,131,3}) \\ &- S(\mathbf{x}\_{2,129}) - S(\mathbf{x}\_{3,130}) - S(\mathbf{x}\_{2,130}) \end{aligned} \tag{9}$$

Here, *x*<sup>129</sup> and *x*<sup>2</sup> appear together in two separate clusters {*x*129,1,2}, {*x*130,129,2}, so their entropy is subtracted. Similar logic follows for (*x*130, *x*3) and (*x*130, *x*2).

We can define the Kikuchi Hamiltonian based on similar intuition. Accordingly, the Hamiltonian for Figure 2B, found using Equation (7), is the following:

$$\begin{aligned} H\_{d=3}(\mathbf{x}) &= H(\mathbf{x}\_{129,1,2}) + H(\mathbf{x}\_{129,130,2}) + H(\mathbf{x}\_{130,2,3}) + H(\mathbf{x}\_{130,131,3}) \\ &- H(\mathbf{x}\_{2,129}) - H(\mathbf{x}\_{3,130}) - H(\mathbf{x}\_{2,130}) \end{aligned} \tag{10}$$

Generally, the approximation accuracy is improved when we consider larger cluster sizes, which is in contrast to mean-field formulation. Explicitly, for *d* = *N*, the Kikuchi approximation for the entropy becomes exact [32]; however, by working with clusters, we cannot define an overall variational distribution, *Q*(**x**), that is consistent with *Qr*(**x***r*) and are unable to obtain an upper bound for *FH* [33].

#### **3. Constructing Cancer Niches using Kikuchi Free Energy**

In this section, we introduce the model used for simulating cancer niches using Kikuchi free energy approximation. For this, we work with the system briefly introduced in Section 2 (Figure 1). Explicitly, this is a 2D Ising model with *N* = 128<sup>2</sup> discrete variables, {*X*1, *X*2, *X*<sup>3</sup> ..., *XN*}, and is arranged in a grid structure. These variables, *Xi*, represent individual cells in a particular tissue population, and each variable can be in one of two states *xi* <sup>∈</sup> {0, 1}. Here, 0 denotes a healthy state, i.e., *<sup>x</sup><sup>o</sup> <sup>i</sup>* , and 1 denotes a cancerous state i.e., *x*<sup>1</sup> *<sup>i</sup>* . We factorise the 2D system as follows:

$$P(\mathbf{x}) = \sum\_{s \in \{0, 1\}} \left( \sum\_{i, j \in N} P(\mathbf{x}\_{i, j'}^s \mathbf{x}\_{i, j+1'}^s \mathbf{x}\_{i-1, j'}^s \mathbf{x}\_{i-1, j+i}^s) \right) \tag{11}$$

where *i*, *j* denote the position on the grid (row, column) and *s* represents the particular variable realisation. This factorisation is a simplified characterisation of sub-populations of cells interacting during cancer niche construction. It is a simple characterisation because it corresponds to interactions of size 4. A different factorization of the system could have been chosen that might require a different higher-order Kikuchi approximation. Technically, even this simple construction leads to non-unique factorisation because of boundary effects. Consequently, to approximate the free energy of this factorised system we need to account for higher-order interactions whilst accounting for overlapping interactions. Specifically, we specify a Kikuchi Hamiltonian to evaluate the energy exchange in these interacting nodes. Practically, this is achieved by using the Kikuchi formulation introduced in for *d* = 3 [14,16] or *B*<sup>3</sup> using [16]'s notation. This is appropriate because *B*<sup>3</sup> with a base cluster as an angle of size 3 gives the same results when using a base cluster of size 4 but with a square [16]:

$$\begin{split} F\_{V} &= \mathcal{U}\_{K}(\boldsymbol{Q}) - S\_{K}(\boldsymbol{Q}) + \lambda\_{1} \left( 1 + \sum\_{s\_{1}, s\_{2}, s\_{3} \in \{0, 1\}} \sum\_{i, j} H(\mathbf{x}\_{i, j}^{s\_{1}}, \mathbf{x}\_{i-1, j-1}^{s\_{2}}, \mathbf{x}\_{i-1, j+i}^{s\_{3}}) \right) \\ &+ 4\lambda\_{2} \left( \begin{array}{c} Q(\mathbf{x}\_{i, i'}^{0}, \mathbf{x}\_{i-1, j-1}^{1}, \mathbf{x}\_{i-1, j+i}^{0}) + Q(\mathbf{x}\_{i, i'}^{1}, \mathbf{x}\_{i-1, j+i}^{1}, \mathbf{x}\_{i-1, j+i}^{0}) + Q(\mathbf{x}\_{i, i'}^{0}, \mathbf{x}\_{i-1, j+i}^{1}, \mathbf{x}\_{i-1, j+i}^{1}) \\ - Q(\mathbf{x}\_{i, i'}^{0}, \mathbf{x}\_{i-1, j-1}^{0}, \mathbf{x}\_{i-1, j+i}^{1}) - Q(\mathbf{x}\_{i, i'}^{1}, \mathbf{x}\_{i-1, j-1}^{0}, \mathbf{x}\_{i-1, j+i}^{0}) - Q(\mathbf{x}\_{i, i'}^{1}, \mathbf{x}\_{i-1, j-1}^{0}, \mathbf{x}\_{i-1, j+i}^{1}) \end{array} \right) \tag{12}$$

$$\mathcal{U}(Q) = \underbrace{\mathcal{U}\_0}\_{=0} + \underbrace{\varepsilon\_{IN}(Q)}\_{Interaction} \tag{13}$$

See Kikuchi and Brush (1976) Table II and IV for graphical representations. Additionally, the Appendix A presents the exact Kikuchi free energy approximation (Equations (A1)–(A3)). Here, *λ*1, *λ*<sup>2</sup> are the Lagrangian multipliers necessary to satisfy the normalisation condition:

$$1 = \sum\_{s\_1 \in \{0, 1\}} Q(\mathbf{x}^{s\_1}) = \sum\_{s\_1, s\_2, s\_3 \in \{0, 1\}} \sum\_{i, j} Q(\mathbf{x}^{s\_1}\_{i, j}, \mathbf{x}^{s\_2}\_{i - 1, j - 1}, \mathbf{x}^{s\_3}\_{i - 1, j + i}) \tag{14}$$

*εIN* is the interaction parameter, and *U*<sup>0</sup> is the activation energy. Following [14,16], we set the activation energy to 0 and quantify the interaction energy using pairwise

interactions (Equation (A2)). Consequently, this particular interaction parameter allows us to accommodate the type of interactions the system evinces. Intuitively, increasing the interaction energy encourages interactions across healthy and cancerous pairs, whilst lower interaction energy encourages interactions within healthy or cancerous pairs [14,18]. Perturbing this interaction energy has a direct impact on the overall free energy of the system. In doing so, we have a way to evaluate the free energy of the system for a given trial. Now, we formalise a setting for how particular realisations evolve (i.e., transition from one state to another) under this system. At each trial, *J* variables may update their current state from either cancer to healthy or healthy to cancer state. This determines the growth rate hyperparameter i.e., *α* = *J N* × 100. The variables transition to the other state (i.e., either healthy or cancerous, given the previous state) determined by the probability *ρ*:

$$\begin{cases} \boldsymbol{v} \sim p(\boldsymbol{v}\_{l}) = \rho^{v\_{l}} (1 - \rho)^{v\_{l}} \\ \boldsymbol{x}\_{j,t} \sim \begin{cases} \boldsymbol{x}\_{j,t} = \boldsymbol{x}\_{\overline{j} \backslash j \in R, t}^{s\_{l}} & \text{if } \boldsymbol{v}\_{l} > \mathbf{n} \mathbf{t}, \\ \boldsymbol{x}\_{j \in N, t}^{s\_{l} \cup s\_{l-1}} & \text{if } \boldsymbol{v}\_{l} \le \mathbf{n} \mathbf{t}, \end{cases} \end{cases} \tag{15}$$

where *vt* denotes an auxiliary Bernoulli random variable indicating whether the variable *xj* is in the same cluster, *<sup>R</sup>*, as *<sup>x</sup><sup>j</sup>* at trial *<sup>t</sup>*, and *nt* is the noise hyperparameter. The variables *<sup>j</sup>* are identified based on the tolerance, *tol*, threshold:

$$x\_{\bar{j}}^{-} \sim \begin{cases} \text{ 1 if } tol > 0.5\\ \text{ 0 if } tol \le 0.5 \end{cases} \tag{16}$$

The tolerance hyperparameter controls the maximum proportion of cancerous states allowed in the system. In our deterministic formulation, these state updates are retained if they minimise the overall Kikuchi free energy of the system. See the Appendix B for the pseudocode (Figure A1).

#### *3.1. Simulations*

Using the above formulation, we simulated three cancer trajectories: local growth, metastasis, and apoptosis. The initialisation of the system is dependent on the simulation, and the specific parameterisations for each simulation are presented in Table 1.

**Table 1.** Parameterisation of the simulations. The interaction parameter regulates the type of cell interaction. The tolerance parameter determines the acceptance level of cancerous cells. The growth hyperparameter controls the number of cell states allowed to switch during a single trial. The noise hyperparameter influences the transition dynamics. These parameters determine the evolution of the system to a steady state according to Equation (14).


#### 3.1.1. Local Growth

First, we modelled growth within the tissue of origin. This speaks to a modification of healthy cells within the original site, i.e., the healthy nodes now become cancerous. Simulating localised cancer cell growth is a core step to elucidating how cancer growth in a healthy organ can be the result of an energy-efficient process in terms of population dynamics although being pathological for the organism as a whole. To simulate this, we specified a high tolerance for cancerous nodes, along with a high energy interaction parameter (Table 1) that constrains the specification of *xj* to variables in the existing cancerous cluster. Here, the positive interaction parameter induces a reduction in overall free energy when cancerous cells interact with other cancerous cells.

For this cancer niche, we initialised the system with a single cancerous cell. Using this, we simulated two distinct patterns: dispersive and localised. The dispersive simulation emulates instances where growth is spread out throughout the tissue population. As a result, the noise hyperparameter was set to 0.25. Conversely, the localised simulation emulates instances where growth is restricted to a small region of the tissue population. For this, we set the noise hyperparameter to 0.

#### 3.1.2. Metastasis

We next modelled metastasis. Understanding metastatic spread is a fundamental priority in cancer research and treatment. With this simulation, we attempted to recreate in silico metastatic progression where the metastatic invasion of a distant organ is enabled by the re-creation of a new niche, without which the cells (although able to migrate) would not survive. Generally, the development of a metastasis occurs in several steps which can develop in a different order over time, including the creation of a premetastatic niche (induced by a distant tumour), mesenchymal transition (EMT) of the original mass (which provides invasive properties), degradation of basement membranes and remodelling of the extracellular matrix (ECM), invasion of the surrounding tissue, angiogenesis, intravasation, arrest of the tumour cells in a capillary bed, extravasation, and the development of macrometastasis [34].

Intuitively, metastasis is the movement of cancer cells from a primary site, (e.g., the original skin tissue of growth) to a secondary location (e.g., lymph nodes). Here, we initialised a small region of the system as the primary site (Figure 3C, grey square), with the remaining grid representing the secondary location. We used a decaying noise hyperparameter to model the transition dynamics alongside an increasing growth rate. This allowed us to emulate the movement of cancerous cells from the primary site to the secondary location, via the bloodstream. Additionally, setting the interaction parameter value to 1.88 allowed us to strike a balance between having both (i) pairwise interactions across cancerous and healthy cells and (ii) pairwise interactions between cancerous and healthy cells within the cluster. This created a system with cancer cells in regions outside the primary site that could support the recreation of a new cancer niche.

#### 3.1.3. Apoptosis

We then modelled apoptosis. This is a controlled mechanism of cell death that intervenes in normal organ repair, growth, and renewal, as well as in inflammation and cancer. It can result from three pathways (extrinsic, intrinsic, and perforine/granzyme). In cancer proliferation, the intrinsic pathway is impaired and cancer cells can escape their death [35]. Here, we aimed to show how the organisms self-preservation mechanisms could be properly activated against the oncogenic cells. Accordingly, we did not expect cancer growth and expected the new niche to be unable to stabilise at a new homeostatic equilibrium.

During this process, the tissue population would normally go through several stages: (i) healthy cells becoming cancerous; (ii) cancerous cells dying; and (iii) remaining healthy cells multiplying and repairing the organ. From our perspective, this simply ensures switching off existing cancerous cells in the population and replacing them with healthy cells.

We were interested in evaluating how the size of the cancer niche impacted apoptosis. Accordingly, we initialised two distinct grids. One had a small cancerous cluster of 41 nodes (or 0.25% of the population), and the other had a larger cluster of 172 nodes (or 1% of the population). For both, we set the tolerance hyperparameter to 10−<sup>5</sup> and the interaction parameter to −1.426 to simulate two instances of apoptosis. Here, having a negative interaction parameter induces a reduction in the overall Kikuchi free energy when more cancer cells interact with healthy cells.

Each parameterisation was simulated across 100 trials with *ρ* = 0.2 and the other parameters were kept consistent.

**Figure 3.** Simulating cancer niches. The figures are pictorial representations of the cancer niches in different trials during simulation. Each row denotes a different simulation: local growth (**A**,**B**), metastasis (**C**), and apoptosis (**D**,**E**). For row C, the grey square represents the original site location.

#### **4. Results**

The results from the simulations, for each cancer niche, are shown in Figures 3 and 4. For the local growth simulation, the construction of the cancer niche was entirely consistent with our expectations. That is, from a single cancerous cell we could observe cancer development. Specifically, by setting a high interaction and tolerance parameter, the cancer was able to survive in healthy tissue. In other words, the overall free energy of the system was minimised by allowing for interactions across cancerous and healthy cells that lead to an oncogenic environment (i.e., the topology of our system) (Figure 4). The free

energy results measured using Kikuchi and mean-field approximation presented similar path trajectories.

**Figure 4.** Free energy as a function of time. The graphic plots the free energies for each cancer niche simulation. Each row denotes a different simulation: local growth (**A**,**B**), metastasis (**C**), and apoptosis (**D**,**E**). For each plot, the Y-axis reports the free energy in natural units [36] and the x-axis as the trial. The first panel reports the Kikuchi free energy (blue dots), interaction energy (red line), and entropy (blue line). The second panel reports the mean-field free energy (blue dots).

The metastasis simulation illustrates (i) the movement of cancerous cells away from the original site and (ii) the ability to sustain the ensuing changes in cell nodes outside the original nodes. (Figure 3). This is a direct result of the interaction parameter value, that allowed for across/between interactions of cancerous and healthy cells. Thus, any changes to the overall system, distal to the original site, were maintained because they minimised the overall Kikuchi free energy of the system. As in the local growth simulation, we see that both Kikuchi free energy and mean-field free energy approximations follow a similar trend (Figure 4); however, the Kikuchi free energy gradient was steeper over time.

For the apoptosis simulation, the system was unable to maintain the existing cancer niche and/or create a new cancer niche. This is reflected by the gradual decline in the overall number of cancerous cells. It is a direct result of reducing the tolerance threshold to 0.00001, where the system now updates the state of cancerous cells to healthy (i.e., a process of renewal and repair). Moreover, the negative interaction parameter value effectuates an overall minimisation of the Kikuchi free energy as the number of interactions between healthy and healthy cells increases. Ultimately, this leads to an apoptotic fate for cancer cells (Figure 4). Interestingly, although we observed a drop in the Kikuchi free energy (as expected), the mean-field free energy for this system increased over time.

#### **5. Discussion**

Organ development and cellular differentiation, de-differentiation, and niche homeostasis arise from a complex interaction between both deterministic and stochastic mechanisms [13]. In these processes, the spatial aspects are as important as the temporal ones, and steady system states are found when the free energy of spatiotemporal dynamics reach minima. In this context, the geometry of the cell's population presents a role that goes beyond the mere effect on mechanics forces [37–40].

The original approach to investigate cellular population dynamics has mainly been deterministic, i.e., a stable genetic code defines cell fate as a programmed hierarchical process. More recently, research has shown how another group of mechanisms based on stochastic self-organisation plays a complementary role in morphogenesis (i.e., organ development) and cell development [41]. Under this perspective, genetic factors would interact with self-organisation mechanisms to balance the growth of an organism and the organisation of cellular assemblies in direct connection with environmental factors. The specific shapes and functions of cells are associated with physical constraints that, although not being genetically encoded, determine the spatial and functional organisation of organs and tissues and further recursively modulate the genetic expression. These physical factors are mechanical, chemical, and geometrical, as well as use-related stressors. Here, the distance or contiguity between cells allows for a spatial-dependent gradient of information such that the regulatory signals created by one cell reach neighbouring cells and regulate the surrounding environment depending on the geography of their position. Therefore, deterministic, stochastic, local, and population factors cooperate during self-organisation, which moves towards a reduction of the overall free energy of the system and stabilisation around the steady states.

Our approach, using Kikuchi free energy, endorses this perspective, namely, that regulatory signals created by cells can affect neighbour cells and influence the surrounding cluster. This induces the creation of new homeostasis in the tissue population (or the body) as a new equilibrium is reached. This stipulates that cancer niche construction is not destructive in physical terms but speaks to the natural evolution of the state function of the system. We have observed changes across our three cancer trajectory simulations. Here, cancer niche construction (or a topographic change to the system) persists because it minimises the overall Kikuchi free energy of the system. This has been shown during metastasis and local growth as nodes change from a healthy to a cancerous state here. Similarly, through apoptosis, a new minimum is reached as the system topology shifts from cancerous to an overall healthy state. The new homeostasis is a consequence of the changes in the interaction and tolerance parameters, that influence the overall Kikuchi free energy estimation and regulate the niche construction.

These parameters allow for changes at both the sub-population (or cluster) level and the overall system. Specifically, under this formulation, the overall system state is regulated by the tolerance level that determines the threshold for having cancerous cells. In our simulations, this meant that a system with high tolerance would permit an increased number of cancer nodes that are conducive to the proliferation and growth of cancer niches. Conversely, low tolerance meant that the system could not maintain the existing cancer niche and/or create a new cancer niche. Additionally, the interaction parameter regulates the types of interactions that would minimise the overall Kikuchi free energy. High interaction parameter values allowed for healthy pairwise interactions for cancer in clusters that are favourable for constructing a new niche. A low interaction parameter value shifted the system state towards either healthy or cancerous by preferring either healthy– healthy or cancer–cancer cell interactions. We postulate that these particular parameters can refer to the biochemical elements in the body which induce cancer formation, growth, and metastasis, as well as the activations of mechanisms of defence.

Our work provides cancer system biology research with a quantitative formulation of evaluating cancer niches that speaks to the recent change of perspective in the field. The role of population dynamics and cancer niches has been proven to be at the core of cancer growth. Our model includes these elements and illustrates their importance in evaluating cancer trajectories. Consequently, this approach has promising future directions for the field of computational biology and can help our understanding of how niches and cells and how possible therapies interact and interfere with each other. Similar mathematical approaches could be considered to test and validate the hypothesis, as well as to predict the possible development of a cancer mass depending on the degree of development and growth within its niche. Although our model is just a reductive simplification of the complex process of carcinogenesis, it demonstrates how Kikuchi free energy is a valuable tool for system biology studies.

#### *5.1. Other Computational Approaches*

The authors of [42] characterised cancer niche concentration as a (stochastic) transition of a healthy system to a distinct oncogenic steady state, e.g., proliferation or apoptosis. They hypothesised that this transition is a direct consequence of the nonlinear dynamic interactions amongst molecular/cellular pathways and modules, e.g., E2F [43], which constitutes an endogenous network. They introduced nonlinear stochastic differential equations (a generalised form of the Langevin equations) in [44–46] to elucidate the specific interactions and nodes that undergird these transitions. See [47] for a review of this approach. This approach is conceptually consistent with a continuous state-space formulation in the current work, albeit introducing stochastic dynamics. Conversely, our approach considers a simplified setting with a discrete-state space formulation with a minimum set of assumptions about the types of nodes and how they interact. Our focus places an endogenous set of agents, as defined by [42], in the setting of cellular population dynamics. This casts the construction of cancer niches in terms of interactions between neighbouring cells and how they influence the cluster using Kikuchi free energy. Interestingly, [47] proposed that free energy can be used to evaluate such noisy systems.

Our approach is complementary to the 2D stochastic cellular automation model proposed in [48] with three states: proliferative, dead ('vacant'), or quiescent. They proposed distinct (deterministic) transitions between each state with three hyperparameters governing the system dynamics (regrowth ability, death rate, and cell cycle arrest). Using Monte Carlo simulations and mean-field phase transition equations, they suggested that the collapse of homeostasis at the multicellular level may be underwritten by non-equilibrium processes; however, their model did not consider long-range intercellular interactions using Kikuchi free energy. Future work could look to incorporate the dynamic transitions introduced in [48] and update rules to model particular cancer niches using Kikuchi free energy.

#### *5.2. Limitations*

There are several limitations with our formulation of cancer niche construction as a consequence of the 2D Ising model used to describe our system. The model restricts the simulations to a closed grid that is unable to interact with the "outside" or be affected by external forces [49]. Nonetheless, this is sufficient for the purposes of understanding how the internal dynamics of a tissue population induce cancer niches. Moreover, our work provides a first step in going beyond a (factorised) mean-field approach to evaluating cancer proliferation and growth [7–10]. A second limitation arises from the deterministic formulation of the transition dynamics. As mentioned above, the self-organisation of cancer niches is a direct consequence of deterministic and stochastic processes that influence the appropriate environment for the growth and maintenance of oncogenic cells. An enactive formulation has been explored in [50] which casts self-organisation as an active inference process where morphogenesis is simply a result of variational free-energy minimization (i.e., of the sort that has been used to explain action and perception in neuroscience [51,52]) and morphogenesis in cellular biology [53].

A key difference between applications of the free energy principle to pattern formation [50] and the treatment in this paper rests upon the nature of the free energy. In applying the free energy principle, the variational free energy pertains to (Bayesian) beliefs parameterised by some (e.g., internal) states about other (e.g., external) states. In contrast, the application of free energy in this paper is directly attributable to the probability of states. This means the free energy principle treats cancer as a process of inference, e.g., a kind of delusion [54], whereas the current treatment treats carcinogenesis as a thermodynamic process (where, under certain conditions, one implies the other); however, future work could look to use Kikuchi free energy under a generalised belief propagation scheme [20] for modelling cancer niches while incorporating (i) external states and (ii) equipping cells with agency by conditioning state transitions on some active states.

Another limitation is a result of the simple generative model, i.e., our discrete random variables can be realised as either cancerous or healthy. Thus, when modelling metastasis, we are unable to model intermediate changes in the cell type (or different realisation of the random variables) as they move from the primary to the metastatic site. Future work could look to expand the model formulation beyond the 2D Ising model to account for these different states as particular cancer niches develop. This would give us a more realistic grounding in the transition dynamics that go beyond being simply healthy to cancerous or cancerous to healthy.

#### **6. Conclusions**

In this work, we illustrate that cancer niche construction is a direct consequence of interactions between clusters of modified cancerous cells using Kikuchi free energy approximation. We show that for certain cancer trajectories, Kikuchi free energy is a more accurate measure of evaluating system topology when compared to a mean-field approximation. Consequently, our work provides proof for the principle of using higherorder free energy approximations that can be more appropriate when evaluating cancer niche construction. Future work should extend the system formulation beyond a 2D Ising construct to evaluate the underlying differences in free energy approximations.

**Author Contributions:** Conceptualization, N.S., L.C. and K.F.; Formal analysis, N.S.; Methodology, N.S., L.C. and K.F.; Visualization, N.S.; Writing—original draft, N.S. and L.C.; Writing—review & editing, N.S., L.C. and K.F. All authors made substantial contributions to conception and design, and writing of the article. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by Medical Research Council (MR/S502522/1, NS), Leverhulme Doctoral Training Programme for the Ecological Study of the Brain (LC), and the Wellcome Trust (Ref: 203147/Z/16/Z and 205103/Z/16/Z, KJF).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

*Exact Formulation for the Kikuchi Free Energy*

We use the Kikuchi free energy formulation introduced in [16]:

$$\begin{split} F\_{V} &= \mathcal{U}\_{K}(\boldsymbol{Q}) - S\_{K}(\boldsymbol{Q}) + \lambda\_{1} \left( 1 + \sum\_{s\_{1}, s\_{2}, s\_{3} \in \{0, 1\}} \sum\_{i, j} H(\mathbf{x}\_{i, j}^{s\_{1}}, \mathbf{x}\_{i-1, j-1}^{s\_{2}}, \mathbf{x}\_{i-1, j+i}^{s\_{3}}) \right) \\ &+ 4\lambda\_{2} \left( \sum\_{i, j; i\_{1}, 2 \colon i\_{i\_{2}}} Q(\mathbf{x}\_{i, j}^{0}, \mathbf{x}\_{i\_{1}, j}^{1}, \mathbf{x}\_{i\_{2}, j}^{0}) + \sum\_{i, j; i\_{1}, 2 \colon i\_{i\_{2}}} Q(\mathbf{x}\_{i, j}^{1}, \mathbf{x}\_{i\_{1}, j\_{1}}^{1}, \mathbf{x}\_{i\_{2}, j\_{2}}^{0}) \right) \\ &+ 4\lambda\_{2} \left( \sum\_{i, j; i\_{1}, 2 \colon i\_{2}}^{0} Q(\mathbf{x}\_{i, j}^{0}, \mathbf{x}\_{i\_{1}, j\_{1}}^{0}, \mathbf{x}\_{i\_{2}, j\_{2}}^{1}) - \sum\_{i, j; i\_{1}, 2 \colon i\_{2}} Q(\mathbf{x}\_{i, j}^{1}, \mathbf{x}\_{i\_{1}, j\_{1}}^{0}, \mathbf{x}\_{i\_{2}, j\_{2}}^{1}) \right) \end{split} \tag{A1}$$

where the following may be represented using Equation (1.3) in [16]:

*U*(*Q*)= *U*<sup>0</sup> - *activation* + *εIN* ⎛ ⎜⎝ ∑ *i*,*j*,*j*1,2,*ii*,2 *Q*(*x*<sup>0</sup> *i*,*j* , *x*<sup>0</sup> *i* 1,*j* 1 , *x*<sup>0</sup> *i* 2,*j* 2 ) + ∑ *i*,*j*,*j*1,2,*ii*,2 *Q*(*x*<sup>0</sup> *i*,*j* , *x*<sup>1</sup> *i* 1,*j* 1 , *x*<sup>0</sup> *i* 2,*j* 2 ) + ∑ *i*,*j*,*j*1,2,*ii*,2 *Q*(*x*<sup>1</sup> *i*,*j* , *x*<sup>0</sup> *i* 1,*j* 1 , *x*<sup>1</sup> *i* 2,*j* 2 ) <sup>−</sup> <sup>∑</sup> *<sup>i</sup>*,*j*,*j*1,2,*ii*,2 *Q*(*x*<sup>1</sup> *i*,*j* , *x*<sup>1</sup> *i* 1,*j* 1 , *x*<sup>1</sup> *i* 2,*j* 2 ) ⎞ ⎟⎠ = *εIN* ⎛ ⎜⎝ ∑ *i*,*j*,*j*1,*i*<sup>1</sup> *Q*(*x*<sup>1</sup> *i*,*j* , *x*<sup>0</sup> *i* 1,*j* 1 ) + ∑ *i*,*j*,*j*1,*i*<sup>1</sup> *Q*(*x*<sup>0</sup> *i*,*j* , *x*<sup>1</sup> *i* 1,*j* 1 )− ∑ *i*,*j*,*j*1,*i*<sup>1</sup> *Q*(*x*<sup>0</sup> *i*,*j* , *x*<sup>0</sup> *i* 1,*j* 1 ) − ∑ *i*,*j*,*j*1,*i*<sup>1</sup> *Q*(*x*<sup>1</sup> *i*,*j* , *x*<sup>1</sup> *i* 1,*j* 1 ) ⎞ ⎟⎠ (A2) *S*(*Q*)= ∑ *s*1,*s*2,*s*3∈{0,1} ∑ *i*,*j S*(*xs*<sup>1</sup> *i*,*j* , *xs*<sup>2</sup> *i*−1,*j*−1 , *xs*<sup>3</sup> *i*−1,*j*+*i* ) − ∑ *s*1,*s*2∈{0,1} ∑ *i*,*j S*(*xs*<sup>1</sup> *i*,*j* , *xs*<sup>2</sup> *i*+1,*j* ) − ∑ *s*1,*s*2∈{0,1} ∑ *i*,*j S*(*xs*<sup>1</sup> *i*,*j* , *xs*<sup>2</sup> *<sup>i</sup>*+1,*j*+1) + ∑*s*1∈{0,1} ∑ *i*,*j S*(*xs*<sup>1</sup> *i*,*j* ) (A3)

#### **Appendix B**

Pseudocode for simulating our system.


**Figure A1.** Pseudocode.

#### **References**

