Unraveling Evolutionary Dynamics: Insights from In Silico Experiments on Selective Mechanisms in Controlled Environments

Unraveling Evolutionary Dynamics: Insights from In Silico Experiments on Selective Mechanisms in Controlled Environments Marco Leddaa, Alessandro Pluchinob, Marco Ragusac aDepartment of Physics and Astronomy, University of Catania, Catania, Italy5 bDepartment of Physics and Astronomy, University of Catania and INFN Catania division, Catania, Italy cDepartment of Biomedical and Biotechnological Sciences—Section of Biology and Genetics, University of Catania, Catania, Italy Abstract10 Driven by the purpose to address the intrinsic difficulties of studying the evolution of biological systems, in this paper, we present a series of in silico experiments aimed at exploring the evolutionary dynamics unfolding in a model of colorectal cancer development, using methods from evolutionary biology. The analysis span from the single gene level, to the phenotypic level. We begin with the standard model commonly used in population genetics to measure natural selection on allele frequencies. We then employ the Price equation, a well-established formalism known for its effectiveness in tracking evolving relationship between entities over time, such as differences between parental and offspring traits. The model results are finally analysed to explain how the selection mechanism works in the context provided by the simulated world. While natural selection is often considered as a statistical phenomenon resulting from changes in genotypic population frequencies, our work suggests that in models with extensive environmental control it is possible to identify the specific elements responsible for exerting selection pressure on the so-called selection units. In this way, we explore the possibility of deciphering the factors that drive and are driven by the selective process. Keywords: Natural selection, Price equation, Colorectal cancer, Evolutionary theory, Computational biology 1. Introduction One of the pivotal goals in modern scientific research is to uncover the dynamical laws or mechanisms governing the behavior of various systems, ranging15 from physical and chemical to biological and social levels. These laws represent the mathematical framework of efficient causation, one of the four Aristotelian Preprint submitted to Journal of L A T EX Templates October 21, 2025 causes accounting for how relatively simple phenomena work.1From this perspective, to describe and predict the behavior of a phenomenon under study involves translation from the qualitative description of the system’s main char-20 acteristics into differential equations. The solution to these equations leads to the deduction of a particular state in the future, or in the past, entirely determined by the premises [1]. In other words, we create a mechanistic and deductive model of reality, which has been an exceptionally powerful method throughout the history of science, discovering and explaining regularities in the25 world around us. However, systems known as complex systems present a particular challenge. They comprise numerous entities from various hierarchical levels, engaging in intricate interactions that defy straightforward explanation through efficient causation. Usually, a linear or nonlinear deductive model is unable to grasp30 the past and future evolution of the aggregate behavior of heterogeneous entities. This complexity is particularly pronounced in biological systems, where entities span different levels, with properties and dynamical relations shifting as one moves from one level to another [2]. Consider, for example, the unclear relationship between structure and function. In fields ranging from molecu-35 lar to evolutionary biology, understanding how a particular function works in conjunction with its corresponding structure does not explain how the same structure can perform new functions in different environments or, conversely, how different structures can perform the same function [3, 4]. Biocomplexity research has highlighted how simple local rules can engen-40 der complex phenomena at higher organizational levels [5]. Therefore, understanding how organisms manipulate matter and exploit energy by deciphering information involves dissecting these entities into their constituent components, studying their unique attributes, and subsequently reconstructing their interrelations following an experimental approach [6, 7]. However, conducting experi-45 ments in fields deeply rooted in history has always been a challenging task [8]. In classical experimental sciences like physics, chemistry, or molecular biology, scientists work with a limited number of variables in a controlled environment, assuming that both endogenous and exogenous perturbations of the system can be integrated into the analysis. Closed environments like labs or individual50 experiments allow for this level of control through the repeated testing of the same process while modifying one or more variables at a time, keeping the other features unchanged. On the other hand, when a historical dimension is inherent to the system, a significant challenge arises: it becomes impossible to intervene multiple times without altering the evolutionary history of the entity in question55 [9]. 1Aristotle articulated the concept of Cause by distinguishing between material, formal, efficient, and final causes. These causes can be understood in relation to the construction of a house: the material cause refers to the physical components used in building the house, the formal cause pertains to the architectural design or blueprint guiding its construction, the efficient cause represents the builder or craftsmen responsible for its creation, and the final cause denotes the ultimate purpose or intention behind constructing the house. 2 Biological systems often span extended time scales with varying rates of change, such as those observed in macroevolution and microevolution (even though the issue is debated [10]). Moreover, once the specific evolutionary history of a system has unfolded, it becomes impossible to restart the process by60 altering one or more variables in isolation to acquire causal insights. In fact, attempting to replay the historical process would yield different outcomes [11]. Nevertheless, a powerful way to comprehend biological systems is to reproduce them, or at least some of their characteristics, with computational methods, where relevant variables are already chosen, and environmental fluctuations are65 more easily manageable [12]. Agent-based models, a specific category of computational models, offer a valuable approach for examining such processes and addressing these issues. These models enable the compression of time scales and the selection of relevant variables, making them more manageable for researchers. Additionally, they facilitate the exploration of spatial dynamics,70 which are essential for understanding evolutionary changes in ecological and population-based contexts [13]. The aim of this paper is multifaceted: first, to utilize agent-based modeling to simulate the evolutionary dynamics of a specific biological system, i.e., the development of colorectal cancer (CRC) [14]; second, to observe how selec-75 tion operates within the system, reasoning about its “causal power” along with unpacking the selective agents; and lastly, to investigate the evolutionary characteristics of the system and how they respond to interventions. The following sections will cover: a short description of the ABM; an analysis of the model using standard population genetics and Price’s equation; a discussion about the80 evolutionary methods in cancer research; issues and open questions resulting from model findings; some thoughts about the agents of natural selection. 2. Methods 2.1. Model’s synthetic description To model the development of CRC within the colonic crypt, we built an85 ABM reproducing the colonic crypt as an unfolded cylinder where 100 layers, of 25 cells each, perform several actions programmed at the individual level [14]. The key elements of the world are two physiological populations, divided into transit-amplifying and differentiated cells, occupying, respectively, the bottom part of the crypt and the top part, in a proportion of roughly 35% and 75%.90 This proportion, as well as the population size, is kept constant by three dynamical factors: the creation of a new bottom layer of cells at each time step; the exit from the simulation of all the cells reaching the top of the crypt; the decrease of β-catenin which starts the differentiation from transit-amplifying to differentiated.95 Each cell has a genome formed by three driver genes (APC,KRAS,TP53 ), which control different functions: APC controls the differentiation of cells and movement; KRAS controls reproduction and movement; TP53 controls reproduction and DNA integrity by means of augmenting the probability of mutation 3 of the other genes when mutated. The three genes can mutate, according to100 the mutation probability, during the mitosis procedure, alone or jointly, either way affecting cellular behavior. Along with genetic features are environmental conditions: an immune-system cell population and a hypoxic threshold. Both conditions are triggered by a population threshold, that is, a given number of cells located in the same spot, which determines the activation of immune-105 system cells or the death of neoplastic cells. All the possible genotypes produce four functional phenotypes upon which we concentrated our analysis, namely: physiological,pre-adenoma,adenoma, and tumoral cells2. The first one, physiological, belongs to cells that have all the alleles in wild-type status or in heterozygosis with no functional changes.110 In the model, genes have two copies represented as lists of binary values ([0/1, 0/1]), denoting the two alleles. These genes are associated with a sum threshold, which is linked to a cellular function such as movement or mitosis time. A value of 0 indicates the wild-type allele, while a value of 1 signifies a mutated allele. For instance, (P2 0APC = 1) if either [1 0] or [0 1]. Cells with phenotypes115 corresponding to physiological,pre-adenoma, and adenoma have a maximum lifespan of 120 hours; tumoral cells have a maximum lifespan of 150 hours. In the following figure, the simulated crypt at different stages of tumor development is shown; see Fig. 1. Runs have been performed for four significant scenarios based on mutation120 probability: µ= 10−9,10−6,10−4,10−2. Three stop conditions: 1 year of model simulation [15]; reversed proportion of transit-amplifying and differentiated cells [16]; population inside of the crypt equals to 10 thousand individuals [16]. 2.1.1. Classic Population Genetics Population genetics is the branch of evolutionary biology dedicated to track-125 ing changes in gene frequencies across successive generations, incorporating fundamental mechanisms like natural selection, random genetic drift, mutation, and migration. These mechanisms exert their influence on biological populations either collectively or independently. Even though there is evidence that neutral theory plays a considerable role in sub-clonal expansion and intratumoral heterogeneity [17, 18], for the purposes of this study, we will focus exclusively on natural selection, as it stands out as the most relevant mechanism when addressing adaptation to novel environmental conditions3[19]. Moreover, its importance as guiding principle of tumor progression has been recognized by several studies during the last decades, see 2In the agent-based model, typology is the variable used to label the four analyzed phenotypes. In other words, we can say that phenotype = typology. 3For our purpose, it is not necessary to model mutation explicitly to observe, for instance, whether or not a mutation-selection equilibrium exists. The reason is that the mutationselection model considers, in a simple two-allele scheme (A1, A2), the reappearance of an allele due to mutation A1←→ A2. However, in our model the wild-type alleles return only from the activity of stem cells, not involving mutation from the neoplastic allele to the wild-type one within the crypt population. 4 [a] [b] [c] Figure 1: Figure [a] shows two local portions of the crypt. The left image depicts the physiological state, where transit-amplifying (blue) and differentiated (green) cells, move upward following the direction of the grey line. The picture on the right depicts a crypt where cells have higher probability of mutation, in fact are present cells with different colors (orange), along with immune-system cells (pink and grey). In [b] we observe, from the left to the right, the four typologies and their different movements and proliferation directions. Solid black lines represent some of the possible movement direction, while dashed black lines are the potential mitosis directions (where the daughter cell is going to appear). Physiological cells have 0°movement, i.e. forward movement; pre-adenoma and adenoma cells have a range spanning from 0°to 90°, both for movement and reproduction; tumoral cells reproduce with a 360° range. In [c] is shown a snapshot with neoplastic cells as to describe the hypoxic effect. The solid-line circle indicates the most superficial cell, and the dashed circle the inner one; when the hypoxic threshold is activated the inner cell is the first one to die and a cascading effect target the others until the value is again below the threshold. for example Fortunato et al. for adaptation [20], Thomas et al. for selection for function [21], or Khong and Restifo presenting natural selection as the escaping mechanism of tumoral cells from immune system [22]. To quantify this variation, the appropriate formalism or mathematical framework is the fundamental model of natural selection [23, 24]. ∆p=pt+1 −pt= ∆p=p(1 −p)(wp−wq) ¯w(1) When this quantity is represented as a function of the wild-type allele fre-130 quency p, it reveals the strength of variation for each allele frequency. If ∆p > 0, natural selection has led to an increase in the frequency of A1in the population. Conversely, if ∆p < 0, the selection process has caused A1to decrease in frequency. Finally, if ∆p= 0, there is no change in the allelic frequency [25]. 5 Although classical population genetics commonly assumes that populations135 exhibit sexual reproduction with random mating, we cannot assume such a condition. Indeed, this feature is absent in our model because the cells reproduce asexually, hence without actually mating with other cells. Regarding the no-recombination assumption, which implies that only complete deletions or point mutations occur, this is also satisfied. The individuals140 (cells) undergo mitosis, producing clones without any genetic exchange with other cells. They are subject to random mutations, but the genetic material remains intact without recombination. The only challenge arises from the assumption of non-overlapping generations, which our model does not fully satisfy. Cells do not die after their first145 mitotic division, so there is always a proportion of cells consisting not only of the contribution of p(or q) to the next generation but also of p(or q) itself. 2.1.2. Price formalism In 1970, George R. Price made a significant contribution to evolutionary biology, and theoretical biology in general, by deriving the equation that bears150 his name [26, 27]. Price’s equation aims to reveal correlations between essential properties (e.g., phenotypes and the number of offspring) at every level of the system under investigation. It serves as a mathematical identity that holds true irrespective of its applicability. Whenever there is a hypothesis aimed at explaining the evolution or change of a biological entity’s trait and its associ-155 ated fitness, this equation represents a valuable tool for uncovering such crucial relationships [28]. Moreover, it has been demonstrated that the most abstract form of this equation is not limited to biological phenomena but extends to processes where information flows from one set to another, with the potential for replication160 errors [29]. Hence, the importance of this formalization of evolutionary dynamics lies in its remarkable generality and simplicity. It is applicable to both asexual and sexual populations, to scenarios involving one or multiple alleles at different loci, and it can account for various factors, such as epistasis, group selection, and multi-level selection [30, 31].165 A wide literature has analyzed the Price equation from both mathematical and philosophical perspectives [32, 24, 33, 34, 29, 35, 36], demonstrating its usefulness [37, 38] and its limitations [39, 40]. A comprehensive derivation of the final form of the equation is provided in the aforementioned literature; here, we aim to show how the equation is constructed and which quantities it explains.170 While Fisher’s model of selection is commonly used in population genetics, the Price equation is better suited to our model since we are not considering sexually reproducing organisms with the assumption of random mating. In our work, we have used this formalism as a framework into which we input data extracted from the model to observe whether changes in phenotypic traits,175 and fluctuations between parents’ and offspring’s phenotypic values, can be explained in terms of selective pressure and their relative importance in fitness. Below, we present one of the most commonly used forms of the equation: 6 ∆¯ ϕ=1 ¯ W[cov(w, ϕ) + E(w¯ δ)] (2) The term ∆¯ ϕon the left side represents the change in the mean phenotypic value. The first term on the right side of the equation, cov(w, ϕ), measures180 how the fitness associated with the phenotype under observation changes over time. This covariance term is quite general, capturing changes caused by both selection and drift. Importantly, it measures a statistical association, reflecting changes directly attributed to differences in the contribution of the ancestral population to the size of the descendant population; that is, it quantifies the185 strength of the selective effect on that trait. The second term, E(w¯ δ), captures the expected differences in the trait values between parents and offspring, scaled by fitness. This term represents changes in the resemblance between ancestors and descendants, which can result from various events both within and outside individuals, such as variations in single190 genes or epigenetic effects (the latter are not included in our model). The index ican refer to any single trait (allele, genotype, phenotype, etc.) or any individual identification belonging to a particular member of the population. Thus, in the present work, if ϕi=xrepresents the trait value of individual lineage iin the parental population, then ϕ′ i=xwill represent the trait value in195 the descendant generation if it remains unchanged, or ϕ′ i=x+nif a mutation has occurred [41]. The linkage between changes in phenotype and subsequent fitness changes is expected in this model, as it reflects how we have programmed the cells to behave; nevertheless, valuable insights can still be gained. To perform an analysis200 of the model using the Price equation, we need to divide the total population into two sub-populations (ancestors and descendants) based on an arbitrarily chosen time interval; see Fig. 2. It is worth noting that this procedure may not capture all possible variations between individuals within and outside of this time interval. One can choose different time intervals, either longer or shorter,205 to examine evolutionary differences between the two sets of populations4. 2.1.3. Genotype-Phenotype mapping In classical population genetics, the fitness of an individual, w, refers to its reproductive success, the contribution of a given phenotype to the next generation, assuming a complete correspondence between genotypic and phenotypic210 traits regarding their contribution to the fitness change [42]. The genotypephenotype mapping can be a complex issue, often represented as a multi-layered graph where various genotypes form the bottom layer, with arrows pointing to one or more phenotypes in the layer above. This complexity arises because one phenotype can be formed by different genotypes, and one genotype can give215 4See the additional material. 7 Figure 2: Here is depicted a possible change in the distributions of ancestor-descendant divided into two sub-populations. The individuals at t1are the ancestors, while the ones at t2are the descendants; different colors and shapes represent different phenotypes. Image built with biorender.com. rise to more than one phenotype, indicating a many-to-many relationship with additional parameters involved [43, 44]. However, in our work the genotype-phenotype map is clearly defined: mutations in both copies of APC yield the pre-adenoma phenotype; adding a mutation in at least one copy of KRAS to the pre-adenoma phenotype results220 in the adenoma phenotype; if both copies of TP53 are mutated, the adenoma phenotype becomes tumoral; see Fig. 3. As a result, the mitosis time only changes in two cases: the mutation of one copy of KRAS and the emergence of the tumoral phenotype. In these cases, the mean mitosis time decreases from 24 hours for physiological and pre-adenoma cells to 12 hours for adenoma cells and225 10 hours for tumoral cells. Following the aforementioned classical definition, we can obtain the proper fitness values by tracking the number of descendants per genotype/phenotype value, hence covering the whole genotype-phenotype space5[45]. Once we have the set of fitness values, we can associate them to the corre-230 spondent frequency plugging them into the recursive equation to track the allele frequency. Here is an example for ∆p= ∆AP C00. ∆APC00 =APC00(1 −APC00)(wapc00 −wapc11 ) ¯w(3) Unlike traditional population genetics, which primarily focuses on measuring changes in allele frequencies within a population, the approach using the Price equation is more general, encompassing everything from changes in allele235 5Again, see additional material for more details. 8 [a] [b] Figure 3: In [a] is depicted the genotype-phenotype map as a table of values. Notice that here are emphasized, with different colors, only the minimum values that trigger the change in cell typologies. However, other combinations present in the chart are functionally different from the physiological one. For instance, in the fourth line the genotype [APC = 00, Kras = 01, TP53 = 00], gives to the bearing cell a reduced mitosis time, because Kras has one mutated copy. Also there are more physiological genotypes other than the [APC = 00, Kras = 00, TP53 = 00]. In the second line we have a genotype, [APC = 00, Kras = 00, TP53 = 01], still in physiological state, since one mutation in TP53 is not enough to trigger different actions. The structure is clearly visible in the plot [b] where is reported the distribution of fitness Win function of the phenotype value ϕ. From the left: the blue rectangle is over the physiological values [0 to 2]; the yellow one is on the pre-adenoma values [2 to 4]; the orange one is on adenoma values [2 to 5]; finally, the red one is on the tumoral values [5 to 6]. frequencies to the characteristics of the entire organism within a population. In our crypt-world model, the phenotypic trait we aim to track across generations is the typology, which is tagged as physiological,pre-adenoma,adenoma, or tumoral; see Fig. 3. The reason behind this choice is the fact that the allelefrequency model is for one gene at a time, hence screening off the possible joint240 contribution of multiple genes to the fitness of the different typologies. In genetic analysis, two types of phenotypic values can be assigned: quantitative and qualitative. A quantitative, or metric, trait is a characteristic of an 9 ations than the non-fixed-time ones. A reasonable explanation for the higher values of ¯ δin non-fixed-time scenarios is that selection maintains traits ensuring higher fitness in fixed-time scenarios. Cells with a beneficial mutation transmit their phenotype to the offspring, and the offspring will proliferate keeping the410 same phenotype, hence reducing the difference between generations after the first one. In non-fixed time scenarios, all phenotypes have the same average mitosis time, thus transmitting differences in movement freedom. However, the difference in movement is not directly linked to a selective advantage, as shown in Fig. 5, hence more variability due to random mutation can occur, increasing415 the difference between offspring and parents. The other interesting phenomenon is the decrease of the term as the mutation probability increases. We can explain this by pointing to the presence of environmental factors targeting neoplastic cells and killing them. Both hypoxic conditions and immune-system cells make the target cellular population size420 drop from the peak, thus resulting in an abrupt change of phenotype values. Nevertheless, what is left is the number of cells of the same lineage with a lower phenotypic value. Suppose one cell with lineage ndivides and, of the two daughter cells, one mutates and the other does not. What happens after a certain period is that the two groups have the same lineage but different phenotypes,425 hence one group will trigger environmental conditions, eventually dying, while the other will survive, lowering ¯ δin the chosen time frame. The change in the mean trait, ∆¯ ϕ, reflects what we derived for the two components of the Price equation, revealing that the overall change, from one generation to another, is greater when a substantial difference in fitness values430 is present. Variance too essentially follows the same trend, showing greater phenotypic heterogeneity as the mutation probability increases [49]. 4. Discussion A fundamental challenge in biological modeling is the attribution of biological meaning to simple mathematical formulations, which are used to isolate435 mechanisms and disentangle causal relations. The mathematical models employed to examine the evolutionary dynamics of the simulated crypt are indeed rather simple; nevertheless, they were able to capture several relevant regularities consistent with evolutionary theory. This “law-seeker” approach does not propose a general model for real tumor evolution, as it lacks several biological440 and ecological features; rather, it aims to offer a general description of the model itself as an abstraction of the real system. The ABM, together with its analysis through evolutionary biology formalisms, has been conceived as a theoretical tool to isolate behaviors and trace back their causes through fine-grained control over both genetic and environmental variables [50, 51].445 However, precisely because of its versatility, the model can help elucidate several theoretical questions, such as ecological interactions or cellular motility as key elements of tumor progression. For instance, it is debated whether predator–prey interaction represents a suitable explanatory framework for the interaction between tumoral cells (prey) and immune system cells (predators)450 16 [52]. Analyzing the dynamics between neoplastic cells and immune-system cells in our model could provide information to confirm or reject the predator–prey model. Invasiveness and metastasis are challenges of paramount importance for cancer therapy and are tightly linked to EMT and the stem-like properties of tu-455 moral cells, which begin to ignore signals from physiological tissues and neighboring cells. Various in silico models have shown how different movement and space-distribution patterns align with experimental findings on metastatic processes [53]. Calibrating the movement degrees and proliferation directions to those of other models would lead to unified frameworks as the basis for more460 sophisticated in silico systems, such as digital twins [54]. Furthermore, as long as certain features of the model can be replicated for the study of different cancers, the model itself can contribute to the theoretical understanding of tumor progression in other tissues and organs. For instance, a statistical association between mutations in TP53 and KRAS has been observed465 in ovarian epithelial cancer [55] and in lung adenocarcinomas [56]. In these cancer types, the activity of both driver genes can be roughly compared to those represented in our model – namely, increased proliferation, altered cell adhesion, and DNA damage leading to a higher mutation rate. Therefore, after adapting the crypt structure to reflect the characteristics of these tissues, it470 may be possible to investigate the roles of these genes in different biological environments. These efforts to build robust, multi-scale models with predictive power naturally complement the possibility of identifying general laws in biology and ecology – defined as statements about invariant properties under given condi-475 tions – which appears plausible if the boundaries of the explanatory system are clearly defined [57]. Thus, once some small-world laws derived from real-world data are established, it would be possible to construct theories about the underlying mechanisms governing their behavior. These mechanisms would, in principle, be able to describe, explain, and predict features of the model, even480 when novel characteristics of the real world are incorporated. For example, if an energy consumption rate linked to a different set of genes were introduced, the same framework developed to study cellular movement should be suitable for explaining the new phenomenon. In our case, that same framework would be the Price equation, accounting for selection and other evolutionary mechanisms.485 Nevertheless, precisely because we are operating within a highly simplified computational model, some issues remain unresolved. Why, for instance, do both setups produce the same trends with only a small difference detected by S? Does that mean that different movements, arising from a given genotype, are the real cause of the difference in the number of mitoses, whereas the fixed average490 mitosis time exerts only a mild influence? Is the fact that wild-type alleles continue to experience positive selection an artifact of the selection model, or do homeostatic mechanisms remain efficient even in such a simplified landscape? These discrepancies could be related to the assumed linearity of the phenomena we are observing [58, 59, 60]. In real systems, such as the colonic crypt,495 a high mutation probability results from a considerable number of factors that 17 strongly deregulate cellular function at all levels, from the genetic to the population scale. These factors cannot simply be summarized and represented by a linear trend, as they influence one another in complex ways, complicating the causal network. Even in an oversimplified model such as the one examined here,500 changes in individual and group-specific traits lead to population effects that, in turn, affect individual traits – specifically, hypoxic threshold cause all the cells in that spot to die even if they are physiological. Otherwise, the phenotypes and their relative fitness would increase in a simple linear manner, explainable by S alone.505 However, as observed by Rice [58], increasing the polynomial degree of the regression equation can indeed improve data fitting, though at the cost of the biological meaningfulness of the values. Figure 7: Three plots showing fit lines for increasing polynomial degrees in red, and the linear fit in yellow. From the top: x2, x3, x4. The linear assumption underpinning the Price equation provides useful, although incomplete, information about the relationship between phenotype and510 18 reproductive success, survival, or any other outcome linked to the perpetuation of a trait across generations. In fact, the first-term coefficient in the linear trend precisely represents the strength of selection on the trait, while the coefficients of higher-order polynomials lose any biological interpretability beyond mere description. In Fig. 7, higher polynomials indicate that something is515 missing from the linear model; however, their coefficients are not related to any biological mechanism (selection) or phenomenon (average fitness). Thus, by balancing simplifying assumptions and biological realism, this modeling approach contributes to the broader challenge for generalizable laws governing complex biological systems.520 4.1. Agents of the selective process Throughout the history of evolutionary biology, the efforts to make evolutionary predictions and to actually predict the evolution of biological phenomena have led to several startling results, highlighting the role of the interaction between genetic and environmental factors [61]. Among these results, much at-525 tention has been paid to the concept of natural selection and its importance in artificial environments [20, 62]. The principle of selection comes into play when it focuses on the transmission of information and the differences in this transmission. In nature, this process occurs as a complex interplay between different organisms struggling for survival, exhibiting both competitive and cooperative530 behavior within and between different species [63]. Moreover, this process is also reproducible in controlled environments, whether in the field, in the laboratory, or in silico [64]. From this perspective, one could legitimately ask: “What are the factors behind this process?” and “Are they decomposable within the hierarchy?”.535 The standard answer to the first question is: a combination of heritable and ecological factors that limit the population’s expansion and influence the individual’s ability to adapt to the environment [65]. If a biological population could reproduce unhindered, it would exhibit exponential, geometric growth that would not be sustainable without limiting factors. Take, for example,540 the bacterium E. coli, which can reproduce every 30 minutes. Given infinite resources, it could produce an astonishing number of copies of itself. Similarly, as Darwin noted, even a single elephant could produce almost 20 million offspring in just 750 years [66]. This exponential growth can also be observed in cell populations, such as545 those that eventually reach a neoplastic phenotype. Early hypotheses about tumor progression assumed that tumor cells proliferate uncontrollably. However, it has been shown that this exponential growth pattern only applies to the initial stages of population growth. In reality, biological populations are limited by several factors, such as the well-known carrying capacity of the system [67]. These550 factors are known as upper-level boundaries or constraints because their effects channel evolutionary or developmental trajectories, constraining the degrees of freedom of lower-level entities. These limits change the dynamic unfolding of the entities’ behavior without altering the nature and range of possible behaviors. 19 An enzyme, for instance, does not change its biochemical capabilities, rather, it555 changes its reaction rate [68]. Let us examine these higher-level interactions within the current model, represented by local population carrying capacity and predator-prey relationships [69]. These factors, along with mutation probability, are responsible for selection pressure on certain phenotypes over others. When we observe population560 dynamics in the presence of environmental conditions, such as cells of the immune system preying on neoplastic cells and a hypoxic threshold causing the death of cells located in the outermost part of the neoplastic mass, we note a tendency toward a homeostatic state characterized by a struggle between cell populations, see (Fig.8 a).565 [a] [b] Figure 8: Crypt simulated at a mutation probability of 10−6, shown (a) with environmental constraints and (b) without environmental constraints. Panel (a) displays multiple neoplastic variants at different crypt locations alongside immune system cells. Panel (b) shows only a single phenotype (orange) due to the absence of environmental factors. Population dynamics plots appear in the top-right inset: green represents differentiated cells; blue indicates transit-amplifying cells; yellow, orange, and red lines correspond to pre-adenoma,adenoma, and tumoral cells, respectively; the pink line represents immune system cells. If, on the other hand, we deactivate the cells of the immune system and increase the hypoxic threshold to thousands of cells, only one population of neoplastic cells explodes in just a few hours (see Fig.8 b). From such a simple control, we can conclude that this dynamic is consistent with the thesis that cells 20 can be considered as individuals of a given species, subject to internal and envi-570 ronmental constraints, and that these environmental elements are consequently the constitutive components of the natural selection process in the model. Given this phenomenology, a conclusion is that: if there are more cellular phenotypes when the immune system cells and the hypoxic threshold are active, therefore is present some kind of causal force acting on individuals and causing575 their actual distribution. Notwithstanding the conceptual problems associated with the term force when applied to natural selection and the evolutionary process [70, 71], there is room for an explanation of natural selection as a causal mechanism rather than just a statistical fact regarding trait distribution. Natural selection is responsible for adaptation to changing conditions; however, the580 causal structure of selective explanation is still a debated issue, with two opposing viewpoints: the statisticalists and the causalists [72]. The statistical view states that selection, fitness, and drift are merely outcomes of the evolutionary process, recognized only at the end of the process. Causalists, on the other hand, argue that these elements are proper causes of the differences in585 organisms’ fitness or changes in the frequencies of organism phenotypes [73]. The standard argument about evolutionary forces has been mainly developed by E. Sober in [74], where he equates natural selection and the other evolutionary forces, to Newtonian mechanics. The case of a null-force is precisely the one described by the Hardy-Weinberg equilibrium, where no forces – either ge-590 netically or ecologically – are acting on the trait, resulting in a null change in the trait distribution. It might be the case of multiple forces acting in opposite directions, resulting in no change at all, due to a balancing effect between those forces, thus connecting for example mutation-selection equilibrium to Newton’s second law. However, if there is a departure from Hardy-Weinberg equilibrium,595 this must be caused by one or more of the known evolutionary mechanisms, detected by the methods developed by evolutionary theory. This mechanistic analysis, while potentially insightful, requires caution regarding whether evolutionary forces genuinely possess the vectorial additivity and deterministic character of physical forces. Nevertheless, this forces frame-600 work provides not a merely metaphorical tool but a genuinely explanatory apparatus, transforming evolutionary biology into a properly causal science where deviations from equilibrium become intelligible as the resultant of identifiable causal factors [75, 76, 77]. Assumed this conceptual framework, we can answer the second question.605 The process of evolution involves the modification of entities’ features across generations through a series of mechanisms, the most relevant of which, in a wide range of events, is natural selection. Natural selection, within this view, is an actual mechanism and a mechanism can be broken down into its component parts to discover the causal role of certain elements in bringing about610 others [78]. The explicit causal elements of the model, in relation to movement and changes in mitosis time, are changes in the lists of genes. A given change in genotype will lead to a corresponding change in the degree of movement of the bearing cell or to a decrease in the average mitosis time. However, when comparing scenarios with fixed versus non-fixed mitosis times, both produce615 21 nearly the same behavioral trends in fitness dynamics. This consistency across scenarios diminishes the relevance of low average mitosis time as a causal element. Although mitosis time and its corresponding fitness changes are triggered by gene mutation, the parallel trends observed suggest that changes in fitness values are predominantly influenced by environmental population factors and620 motility factors, rather than by mitosis time itself - probably, even in scenarios where mitosis time is predetermined. If we look at the plots in (Fig. 5), we can observe that in the absence of a specific decrease in the mean mitosis time when the genes are mutated, the fitness values still increase up to the maximum value at ϕ= 6. This trend is625 almost the same as in the simulations with predefined changes in mitosis time. Therefore, changes in movement, which can be considered as different degrees of freedom for the cells, or equivalent to the absence of individual constraints, are always relevant for the number of offspring produced. In this sense, we can argue that population-level constraints are as relevant as changes in mitosis time630 as causes of a selective process. A major challenge also arises in the attempt to define the entities subject to selection, a problem that has puzzled biologists since Darwin’s time [32]. In the hierarchical view of biological systems, entities are not isolated; they are nested in collectives that have structures and functions distinct from those of their635 constituent units [79]. Cancer is usually viewed as a multilevel phenomenon, ranging from the structures of the cell nuclei to the tumor microenvironment, for which there is still no clear spatial and relational definition [80]. Consequently, natural selection occurs at multiple levels of organization in cancer, not only at cellular and organismal levels, but also at the level of genes (includ-640 ing extrachromosomal DNA and retrotransposons), mitochondria, collections of cells such as epithelial proliferative units, and potentially even metastases and organisms [62]. In this view, the roles of transit-amplifying and differentiated cells are to keep the crypt population in balance, secrete mucus, and absorb nutrients from645 substances flowing down the gut. The function of the entire crypt, and of the colon tissue as a whole, is to ensure various processes that result from the coordinated behavior of multiple cell types [81] and even bacterial populations [82]. Therefore, the structure and functions of both the crypt and the gut, which are of great importance in the hierarchy, act as selective agents for the behavior650 of individual cells. At the population level, these behaviors, in turn, become selective elements for higher-level structures and functions. These latter aspects are the macroscopic result of individual behavior at the microand mesoscopic levels. Thus, when cells acquire somatic mutations that alter their behavior, they are referred to as cheater or renegade cells because the selective process no655 longer benefits the entire tissue but rather the individual unhindered cell. For these neoplastic variants, the collective behavior no longer serves as a fitness enhancer [83]. 22 5. Conclusions In this study, we performed several simulations with the aim of gaining valu-660 able insights on the existence of natural selection mechanisms and how these mechanisms drive the behavior of both individual cells and cell populations. Interesting conclusions can be drawn from the results obtained when predetermined mitosis time was not used, i.e. different degrees of freedom in movement lead to different fitness values. It is biologically reasonable to assume that for665 physiological behavior at the individual and population level, tissue architecture (movement allowed in one direction only) is a limit for the individual cellular activity. Hence, when individual mutations deregulate local sub-populations, macro-level boundaries are altered or are turned ineffective in the limitation of individual behavior (allowing movement in all directions), thus leading to670 abnormal and dangerous phenomena such as uncontrolled mitosis. A further conclusion can be drawn from the results indicating that environmental conditions play a relevant role in the interplay of different neoplastic phenotypes. These interactions are not present in the simulations of the crypt, in which the cells of the immune system and the hypoxic thresholds are inactive. From that675 we can reasonably deduce that, when considering an evolutionary perspective of tumor progression, we must consider environmental factors not only as background conditions, but as true causal (a component of the selective mechanism) factors alike the activation of oncogenes and the silencing of tumor suppressor genes by mutational events.680 We believe that these results place the model in the third category of Levins’s classification [84], which is to sacrifice precision for realism and generality. A certain degree of realism is implicit in the model because, although it is based on simplified assumptions, it relies on experimental and clinical data. Therefore, as long as the importance of driver genes in tumor development is confirmed685 by the majority of clinical samples, a small part of reality is captured by the simulated world. Generality refers to the importance of local behavior in generating emergent organizational properties that counteract shaping the evolution of the whole system. In the simulated world, the local behavior is characterized by the movement and the reproduction direction of the agents, which are690 almost universally shared properties among different living beings [85]. These local features in turn triggers population thresholds which cause either death or opportunity for the individuals belonging to different cellular populations [86, 87]. Our future research efforts aim to formally define the causal evolutionary695 structure underlying the model and to clearly dissect the selective agents to gain insights into their weight in the progression of CRC. Moreover, these results will also be useful to shed light on the concept of top-down causal relationships or higher level constraints and how they can be formalized using the language of causal inference [88, 60].700 23 6. Appendix Datasets and computational model are available at the following links: Google documents (necessary for the size of the data), https://drive. google.com/drive/folders/1iLLbqHcUKkQKO_qRvBZ8I92IGspZdF5_?usp= sharing705 Github, https://github.com/MLedda7/CRC-data-models References [1] S. A. Kauffman, Prolegomenon to patterns in evolution, Biosystems 123 (2014) 3–8. [2] A. Wagner, Causality in complex systems, Biology and Philosophy 14710 (1999) 83–101. [3] G. F. Ellis, Efficient, formal, material, and final causes in biology and technology, Entropy 25 (9) (2023) 1301. [4] M. Herman, B. Aiello, J. DeLong, H. Garcia-Ruiz, A. Gonz´alez, W. Hwang, C. McBeth, E. Stojkovi´c, M. Trakselis, N. Yakoby, A unifying framework for715 understanding biological structures and functions across levels of biological organization, Integrative and Comparative Biology 61 (6) (2021) 2038– 2047. [5] P. Hogeweg, Multilevel cellular automata as a tool for studying bioinformatic processes, in: Simulating complex systems by cellular automata,720 Springer, 2010, pp. 19–28. [6] F. J. Bruggeman, H. V. Westerhoff, The nature of systems biology, TRENDS in Microbiology 15 (1) (2007) 45–50. [7] H. Kitano, Computational systems biology, Nature 420 (6912) (2002) 206– 210.725 [8] R. E. Lenski, M. Travisano, Dynamics of adaptation and diversification: a 10,000-generation experiment with bacterial populations., Proceedings of the National Academy of Sciences 91 (15) (1994) 6808–6814. [9] B. Batut, D. P. Parsons, S. Fischer, G. Beslon, C. Knibbe, In silico experimental evolution: a tool to test evolutionary scenarios, BMC bioinformatics730 14 (Suppl 15) (2013) S11. [10] A. M. Simons, The continuity of microevolution and macroevolution, Journal of Evolutionary Biology 15 (5) (2002) 688–701. [11] S. J. Gould, The structure of evolutionary theory, Harvard university press, 2002.735 24 [12] D. C. Krakauer, J. P. Collins, D. Erwin, J. C. Flack, W. Fontana, M. D. Laubichler, S. J. Prohaska, G. B. West, P. F. Stadler, The challenges and scope of theoretical biology, Journal of theoretical biology 276 (1) (2011) 269–276. [13] N. Takeuchi, P. Hogeweg, Multilevel selection in models of prebiotic evo-740 lution ii: a direct comparison of compartmentalization and spatial selforganization, PLoS Comput. Biol. 5 (10) (2009) e1000542. [14] M. Ledda, A. Pluchino, M. Ragusa, Exploring the role of genetic and environmental features in colorectal cancer development: An agent-based approach, Entropy 26 (11) (2024) 923.745 [15] T. Ingham-Dempster, B. Corfe, D. Walker, A cellular based model of the colon crypt suggests novel effects for apc phenotype in colorectal carcinogenesis, Journal of computational science 24 (2018) 125–131. [16] K. Mamis, R. Zhang, I. Bozic, Stochastic model for cell population dynamics quantifies homeostasis in colonic crypts and its disruption in early750 tumorigenesis, Proceedings of the Royal Society B 290 (2009) (2023) 20231020. [17] M. J. Williams, B. Werner, C. P. Barnes, T. A. Graham, A. Sottoriva, Identification of neutral tumor evolution across cancer types, Nature genetics 48 (3) (2016) 238–244.755 [18] A. Niida, K. Mimori, T. Shibata, S. Miyano, Modeling colorectal cancer evolution, Journal of Human Genetics 66 (9) (2021) 869–878. [19] B. Ujvari, B. Roche, F. Thomas, Ecology and evolution of cancer, Academic Press, 2017. [20] A. Fortunato, A. Boddy, D. Mallo, A. Aktipis, C. C. Maley, J. W. Pepper,760 Natural selection in cancer biology: from molecular snowflakes to trait hallmarks, Cold Spring Harbor perspectives in medicine 7 (2) (2017) a029652. [21] F. Thomas, J. DeGregori, A. Marusyk, A. M. Dujon, B. Ujvari, J.-P. Capp, R. Gatenby, A. M. Nedelcu, A new perspective on tumor progression: Evolution via selection for function, Evolution, Medicine, and Public Health765 12 (1) (2024) 172–177. [22] H. T. Khong, N. P. Restifo, Natural selection of tumor variants in the generation of “tumor escape” phenotypes, Nature immunology 3 (11) (2002) 999–1005. [23] C. R. Linnen, H. E. Hoekstra, Measuring natural selection on genotypes and770 phenotypes in the wild, in: Cold Spring Harbor Symposia on Quantitative Biology, Vol. 74, Cold Spring Harbor Laboratory Press, 2009, pp. 155–168. 25