Network Inference from Cancer and Epithelial Cell Co-Culture Reveal Microenvironmental Configurations That Drives Population Dynamics

1. Introduction

In cell co-culture experiments, after eight days, keratinocytes surround the melanoma clusters, limiting the growth of the melanoma cells. The keratinocytes are denser near the melanoma clusters and less dense further away [1]. It’s interesting that when keratinocytes are cultured with melanoma, they become less differentiated and connect with the melanoma instead of other keratinocytes nearby [2,3]. These observations may explain biologically the different levels of tolerance to high density between melanoma and keratinocytes [1].

As in concepts about continuous/discrete dynamic system modeling, considering the time interval, this spatial characteristic observed at the end of the experiment corresponds to an equilibrium state or attractor of the population dynamics. Conversely, the initial distribution of cells on the cell culture plate is the initial condition of the basin of attraction. Thus, the system shift from the initial condition to its attractor through growth dynamics such proliferation, migration, survival, and cell death, causing changes in cell density and in the spatial pattern of the population.

Herein, among system biology approaches to modeling population dynamics is the Agent-Based Model (ABM) and Boolean Networks (BN). ABM is a bottom-up simulation technique where agents interact with each other [4]. Among the ABM classes, one can find the Cellular Automata (CA), which is a collection of cells on a grid that evolves through a number of discrete time steps, obeying rules based on the states of neighboring cells, such that a spatial patterns emerge [5,6]. Firstly, the CA was one-dimensional, where the rules of cell state transition were dependent on cells alignment in the neighborhood, but in the first two-dimensional CA, the Conway Game of Life, also known as a life-like CA, the rules of state transition of the cells were dependent of neighborhood cell count and described as inequalities [7]. The Square-Lattice Cellular Automata (SLCA) in turn, a two dimensional CA, has been successfully used to simulate avascular and vascular tumor growth [12,13,14,15,16,17,18,19], invasion [20,21,22,23], tumor interactions with environmental cues [24,25,26,27,28,29,30], stromal composition [30,31,32,33,34] and tissue architecture [35,36,37,38].

On the other hand, BN are models wide used in gene expression or signaling pathway studies. They consist of a collection of variables (the nodes), with binary states determined by other variables in the network through logical functions [8,9]. There is already efforts to link BN and population dynamic such as microbiome gut genus abundance through BN inference [10,11], and they reveal ecological interaction among bacterias. Similarly, BN has been applied in many fields, including quantitative System Pharmacology [39], competition in microbiome [10], stability of apoptosis [40,41], epithelial-mesenchymal transition [42] and target therapy [43].

CA is also considered a network because each neighboring cell is wired to a focal cell and vice versa, and the rules determined by the local dynamic of all N cells in the grid contributes to the evolution of CA [46]. Based on this, some authors inverted the forward mappings of CA local dynamics to calculate the backward trajectory, thereby unveiling the network of CA dynamics also known as global dynamics [44,45,46,47,48,49]. For this reason, CA are generalizations of BN one another in terms of basins of attraction and are considered instances of discrete dynamical systems composed of parallel-acting components [46]. Here, we argue that if one can find a network underlying an ABM dynamic such CA and reveal its basins of attraction in a BN manner, one can estimate rules of an ABM evolution based in neighborhood cells state through the BN state transitions within the basin of attraction.

We argue the above because the cell density is a fundamental dimensional unit in a cell culture plate, and one can qualitatively characterize the cell density state as "high" or "low", like the confluence itself. So BN inference from co-culture growth curve data would provide the cell density state that changes over time. By comparing these cells density states within state transitions and along a basin of attraction, we can derive an inequality in relation to cell neighborhood count, such as life-like CA, that can describe any regulation among cells because the growth curve emerges from local regulations amid the Tumor Microenvironment (TME) [50,51,52,53]. We aim (1) to infer a BN from the growth curve time series obtained from in vitro co-culture, (2) to estimate the rule ABM from the IBN basin of attraction by comparing the density states within state transitions and reaching a consensus among state transitions about density states, and (3) to implement an ABM to simulate spatial dynamics.

We observed that the binarization of the time series points to a high population density after the fourth day. The IBN functions allowed us to choose a regulator for each node that is supported by the literature, and the rules estimated from the IBN basin of attraction point to conditions for cells migration and proliferation in agreement with experimental data [1,3]. The rules aforementioned do not embrace known conditions to cells survival and death, therefore we ended up proposing intrinsic conditions to cells survive and die. The migration of keratinocytes increases smoothly until the fourth day of the simulation, making them denser near melanoma clusters. In the mean time, the survival probability increases substantially, constraining melanoma proliferation. Therefore, the simulation suggests that the attractor state of cell co-culture population dynamics could be induced by an increase in keratinocyte migration and survival.

2. Material and Methods

In this section, we present the methodological approach used for data inference and modeling. We shall highlight that the modeling process is contingent upon the outcomes derived from the inference. Figure 1 illustrates the workflow for data analysis and modeling. The source code link can be found at S4 Source code.

2.1. Experimental Data

We use previously published cell proliferation data measured in co-culture [1]. The data are the density ($cells/m{m}^2$) of normal human keratinocytes (HaCaT), metastatic melanoma (SK-MEL-147), at co-culture with the initial proportion of 1:10 (SK-MEL-147:HaCaT) [54] measured along 8 days in triplicate. To quantify the cell density necessary to achieve the carrying capacity, we introduced the variable S to express the difference between the carrying capacity (K) and the population density at each specific time point. This additional variable is justified because the space has a significant impact on the population dynamics. Hence, the more cells need to reach carrying capacity, less space is available, and the fewer cells need to reach carrying capacity, more space is available. Consequently, this variable serves as an initial arithmetic approximation of the space availability. Throughout this paper T will represent melanoma, H keratinocyte, and S space availability.

2.2. Boolean Network Inference

BN inference involves two steps: 1) binarization, which transforms cell counts per day in binary discrete data, and 2) network learning. We evaluated three different binarization methods: (1) k-means, (2) iterative k-means [10,55] and (3) BASC [56]. The binarization outcome chosen to network learning shall be the one which point highest density after confluence onset in agreement with Morais et al [1]. We used BoolNet 2.1.8 from R 4.2 [57] to implement the binarization by BASC and k-means, and Python 3.10.2 to implement the binarization by iterative k-means, adapted from Steinway et al [10]. In the context of this work, the number 0 denotes low density, and the number 1 denotes high density. We performed BN learning using the BestFit algorithm [58] through the BoolNet package 2.1.8 from R 4.2 [57].

2.3. Estimating the Agent-Based Models Rules from BN Basins through Consensus of State Transition Comparison

Here propose a simple method to estimate ABM rules from BN basin of attraction. The proposed approach involves two steps. First, we compare the density states among variables within state transitions belonging to a given basin of attraction and generate some sub-logic rules in term of inequalities. Second, we obtain a consensus about sub-logic rules among all state transitions belonging to the same basin of attraction. For instance, while analyzing six state transitions, it is observed that a specific variable A = 1 and variable B = 0 ($A\ge B$) in three of them, while A = 0 and B = 1 ($A\le B$) in one of them, and A = 0 and B = 0 ($A=B$) in two of them. The most observed inequalities are $A\ge B$ and $A=B$, and we define the ABM rule as $A\ge B$. We determine absence of regulation between variables when >, = and <, have the same frequencies such in case.

2.4. Implementation of Agent Based Model

Among ABM classes we’ve implemented a SLCA which was built based on the estimation of growth dynamic rule from the BN basin (Figure 1). The model consists of a tissue scale, resulting in a single discrete variable, namely the cells. The model encompasses a dynamic process involving cell proliferation, death, migration, and survival. This process runs on a two-dimensional grid ($n\times n$ with n = 200) with Moore neighborhood. The criteria for determining grid size was based on the spatial model developed by Morais et al [1]. Each vertex can accommodate only one tumor cell. During each iteration, all cells undergo the actions outlined in Figure 5. The cells were randomly distributed over the grid in a ratio of 1 T to 10 H according to [1]. We assume a time-step of 12 hours resulting in a total of 16 iteration, corresponding to 8 days-culture. Section 3.3 describes further detail about the state transition rules.

2.5. Spatial Configuration, Parameter Estimation, and Population Dynamics Process Probabilities

We applied a convolutional Gaussian filter of size $k=5\times 5$ and standard deviation ${\sigma }^2$ to estimate the relative local density of H cells (arbitrary units) in the simulation domain. To calculate the growth curve of both cell lines we fitted the logistic model (Equation (1)) using the curve fit function of the optimize module of the Scipy library (version 1.7.1) in Python (version 3.10.2)

$$N\left(t\right)=\frac{K}{1-e^{-\rho (t-\tau )}},$$

(1)

where K denotes carrying capacity, $\rho$ denotes proliferation rates and $\tau$ denotes time for half asymptote.

Throughout the simulations, we calculate the probability of the cells going through those tumor growth dynamic processes described in Section 2.4. Denote $i=1,...,n$ as an element of the list $\omega$ = (Proliferation, Death, Migration, Survive) and f the frequency of occurrence of one of the processes in the list $\omega$, thus

$$Pr\left({\omega }_i\right)=\frac{f_{{\omega }_i}}{{\sum }_{i=1}^n{f}_{{\omega }_i}},$$

(2)

for both T and H.

3. Results

3.1. The Iterative k-Means Binarization Method Demonstrated That Both Populations Are in High Density after 4 Days

As described in Section 2.2, the time series was binarized, and then we learned a BN. The binarization by iterative k-means showed a high density of H and T and a low availability of S, after four days (Figure 1 and in Figure S1).

Figure 2. The heatmap represents the binarization of co-culture time series by iterative k-means. The yellow cells represent 1, which means high density, and the blue cells represent 0, which means low density. Each row represents the binarized trajectories of H, T, and S, and each column represents the state transition of the network over eight days.

Figure 2. The heatmap represents the binarization of co-culture time series by iterative k-means. The yellow cells represent 1, which means high density, and the blue cells represent 0, which means low density. Each row represents the binarized trajectories of H, T, and S, and each column represents the state transition of the network over eight days.

Table 1 shows the possible Boolean functions for each variable. The possible function for each variable is highlighted in yellow with its explanation. The intention of showing two different biologically reasonable networks is to ascertain possible plasticity in the estimated rules for ABM as the BN basin topologies change once a BN has $2^n$ possible state transitions, here n=3 (H, T, and S).

3.2. The Network Dynamics Exhibit Two Singleton Attractors

After we binarized the time series and applied a method of network learning, we generated the state graph, and observed two singletons (fix points) attractors in both IBN 1 (Figure 3a) and 2 (Figure 3b). Before evaluating the attractors and proceeding with the estimation of CA dynamics rules, we labeled each transition state as $B_1,B_2,B_3,B_4,B_5,B_6,A_1,A_2$. These IDs belongs to the same state transition in both IBN 1 and 2. The attractor $A_1$ shows T and H in low density, while S is in high density, which ultimately means high space availability or the population is scattered. The attractor $A_2$ shows H and T in high density and S in low density, which ultimately means low space availability or that the population is in high density. The IBN 1 has no predecessor state transition from initial condition among basins of attraction, while the IBN 2 present an small sequence of state transition not.

3.3. Boolean Network Basin Provide Growth Dynamic Rules for Proliferation and Migration

To estimate the growth dynamic rules to ABM from BN basin of attraction, we analyze the state transitions in terms of density state comparison. We ignore the state transitions that do not have a reasonable biological meaning considering a cell culture plate. $B_3$ presents a low cell density and low space availability, and $B_6$ presents a high cell density and high space availability. When we compare the density states between T and S among state transitions $B_1$, $B_2$, and $A_1$, $T<S$ is in two-thirds of these state transitions (orange-dash rectangle at Figure 4). Additionally, $T<H$ is in 2/3 of the aforementioned state transitions (orange-dash rectangle at Figure 4). Thus, the consensus of the density state comparisons among these state transitions shows two inequalities: $T\le S$ and $T\le H$ (Figure 4). In contrast, upon comparing the density states between T and S among transition states $B_4$, $B_5$, and $A_2$, $T\ge S$ is in two-thirds of these state transitions (purple-dash rectangle at Figure 4). Furthermore, $T\ge H$ is in 2/3 of those state transitions (purple-dash rectangle at Figure 4), so by consensus among these state transitions, we obtained two new inequalities: $T\ge S$ and $T\ge H$ (Figure 4).

The sub-logic rules $T\le S$ and $T\le H$ represent a condition in which there is a low amount of T and a high amount of H or S, which do not favor an increase in density, and according to [3], we determine that it favors migration. Thus, we set Equation (3) where $M_{T,H}$ denotes migration.

$$M_{T,H}=(T\le S)\vee (T\le H)$$

(3)

The sub-logic rule $T\ge S$ and $T\ge H$ represent a condition where a higher density of T and a lower density H or S, favors density increase. Thus, $P_T$ and $P_H$ denotes proliferation of T and H in Equations (4) and (5), respectively. Due to the lack of data about the resources present in the cell culture media, we decided to augment Equations (4) and (5) with the probability of instantaneous proliferation $1-\frac{N}{K}$[60], since the logistic model fitted to the growth curves present in Morais et al [1]. If the probability of instantaneous proliferation is satisfied ($1-\frac{N}{K}\ge r$, where r is a random number) then the cells can proliferate. Another important condition to simulate cell dynamics is cell death, represented by the death term ${\overline{D}}_T$ in Equation (4) and (5), and the survive term ${\overline{SV}}_H$ that blocks H cell proliferation, as explained in Section 3.4.

$$P_H=\left(1-\frac{N}{K}\vee ((T\ge S)\vee (T\ge H))\right)\wedge {\overline{SV}}_H$$

(4)

$$P_T=\left(1-\frac{N}{K}\vee ((T\ge S)\vee (T\ge H))\right)\wedge {\overline{D}}_T$$

(5)

3.4. Data from Cell Lines and Growth Curves Allow Us to Propose an Intrinsic Rule to Describe Survival and Death of Each Agent

Here, we only could estimate a rule based in neighbors count that are compatible with migration and proliferation. Hence, we assigned an intrinsic rule to the agents survives or die. So in order to describe the survival and death of cells, we turn to Hayflick’s limit [61,62,63] concept to H and cancer stem cell hypothesis [64,65,66,67,68,69] to T because both can embrace the so-called parameter proliferation capacity (${\rho }_{max}$). Thus, with each cell division, the proliferation capacity of T and H will decrease by one unit, ${\rho }_{max}={\rho }_{max}-1$. It is well documented that H undergoes a limited number of duplications in the cell culture plate and does not always reach 100% confluence; thus, once H reaches ${\rho }_{max}=0$, they survive ($S{V}_H$), as proposed in Equation (6). We assume that when T reaches ${\rho }_{max}=0$ and in circumstances where there is no space ($D_T$), the crowding situations might yield cell de-adherence and death [70], as in Equation (7).

$$S{V}_H={\rho }_{max}=0$$

(6)

$$D_T=\left({\rho }_{max}=0\right)\wedge \left(S=0\right)$$

(7)

Once H does not migrate, proliferate, or survive, they die, as in Equation (8). Alternately, once T does not migrate, proliferate, and die, they survive, Equation (9).

$$D_H={\overline{M}}_H\wedge {\overline{P}}_H\wedge {\overline{SV}}_H$$

(8)

$$S{V}_T={\overline{M}}_T\wedge {\overline{D}}_T\wedge {\overline{P}}_T$$

(9)

Figure 5. Flowchart of the ABM implementation. First, we set the initial density of cells and a random initial distribution of them over grid. At each iteration, the algorithm updates the amount of T, H, and empty vertex (S) in all focal cell neighborhoods. The conditions described in Equations (3) to (9) are updated in each cell. When there is at least one empty vertex, the cells proliferate or migrate, but when the proliferation capacity is exhausted, the cells survive or die.

Figure 5. Flowchart of the ABM implementation. First, we set the initial density of cells and a random initial distribution of them over grid. At each iteration, the algorithm updates the amount of T, H, and empty vertex (S) in all focal cell neighborhoods. The conditions described in Equations (3) to (9) are updated in each cell. When there is at least one empty vertex, the cells proliferate or migrate, but when the proliferation capacity is exhausted, the cells survive or die.

3.5. The Parameter Estimation Are in Agreement with Experimental Data

The values of cells density and the frequency of proliferating, migrating, surviving, and dying cell in the overall simulation were used to evaluate if the growth curve fit to the logistic model (Equation (1)). The growth curves fit the logistic model (Figure 6g) and the estimated values for the parameter K, $\rho$, and $\tau$ are presented at Table 2. Considering each iteration as a $\Delta t=12$ hours favored the estimation of ${\rho }_T\approx {\rho }_H$ and ${\tau }_T\approx {\tau }_H$. We shall recall that the results in Figure 6 concern ${\rho }_{ma{x}_T}=7$ and ${\rho }_{ma{x}_H}=6$. According to Morais et al [1] the values estimated to ${\rho }_T$ and ${\rho }_H$ are equal as well those estimated to ${\tau }_T$ and ${\tau }_H$.

3.6. The Density of H Indeed Is Higher Close to T Cluster

We implemented the spatial model using ABM (Section 2.4, Section 3.3 and Section 3.4) and, starting from a given an initial distribution of cells. on the fourth day we can already see some T clusters surrounded by H, so this spatial pattern is consolidated on the eighth day (Figure 6a–c and animation in S3 Video). The heat map shows H density in the neighborhood of each focal cell (Figure 6f) points to higher yellow intensity around T clusters. This can be explained by the smooth increase in H migration probability from the beginning of the simulation (Figure 6(h)) as well as the increase in T survival probability from day 4 onward (Figure 6g). Therefore, as T clusters consolidate, more of the H neighborhood lacks T; hence the probability of migration decrease mostly because the probability of H survive increase substantially after four days (Figure 6(h)).

4. Discussion

In this study, we used the time series of the growth curve in the co-culture of T and H to infer a BN. Through BN inference, we were able to determine the boolean functions that led to network dynamics. These dynamics were shown in a state graph, where we arranged the state transitions according to the basin of attraction and estimated the ABM rules. These estimated rules were used in the implementation of an ABM to simulate spatial population dynamic. The binarization of the growth curve data reveals that the population attained a high density after 4 days. This observation aligns with the findings of Morais et al [1], who reported that the population achieves confluence after the same 4-day period. The proposed methodology captured the relationship between the confluence observed on the cell culture plate and the corresponding cell density. A boolean network of n nodes exhibits $2^n$ state transitions; hence, our time series, consisting of a sequence of 8 data points, can fully define the system. The occurrence of non-null error rates in the inferred boolean network arises due to the similarity of the H and T binarized trajectories, also known an identifiability, nevertheless, this issue can be solved by prior biological knowledge or by pruning the network nodes according to the context [39]. The BestFit BN learning paradigm searches for possible boolean functions that describe the state of a given variable (node) at time $t+1$ as a function of the network state (the states of all nodes) at time t. In the heatmap (Figure 2), S is ON when both T and H are OFF, and vice versa, which means that S inhibits both T and H, which in turn inhibits S, and finally both T and H can "activate" each other. The IBN 1 in Table 1 and Figure 3a shows a canalizing node and a initial condition with no predecessor state transitons throughout basin of attraction. The IBN 2 in table Table 1 and Figure 3b, on the other hand, shows a small number of states before reaching the attractor. However, both networks (Figure 3a,b) have the same attractor. Hence, different BNs can lead to the same attractors. This indicates that the estimation of ABM rules was not affected by changes in the topology of the state graph, suggesting a system plastic capacity, as in the genetic networks [71]. Indeed melanoma affects keratinocytes phenotypes [2,3], our approach only captured an attractor state concerning local density of both melanoma and keratinocyte. By comparing the density states within state transitions $B_4$, $B_5$, and $A_2$, we could estimate $(T\ge S)\vee (T\ge H)$, suggesting that T will only proliferate where there is more T, since experiments with model organisms like zebrafish have shown that melanoma tends to cluster [72]. Likewise, we could estimate $(T\le S)\vee (T\le H)$ based on how the density states within states transition $B_1$, $B_2$, and $A_1$, meaning that either H or T will move away [1,3]. It is worth highlighting that these Boolean expressions can be concurrently satisfied in a dynamic context, so the potential for cell proliferation and migration is inherently incorporated, which reduces any challenges in modeling associated with the go-or-grow phenomenon [73,74,75,76,77]. The spatial dynamic shows H more denser nearby T clusters than close to themselves (Figure 6a,b); moreover, the values estimated for ${\rho }_{T,H}$ and ${\tau }_{T,H}$ are very similar. The value estimated to ${\tau }_{T,H}$ is initial cell density-dependent and it is in agreement with [1]. The growth curve certainly did not emerge from the dynamics itself, but it is likely hardwired to the instantaneous probability term present in Equation (5) and (6) that was employed to represent the resources in the culture medium. Normal cells like H undergo a limited number of duplications in the cell culture plate and do not always reach 100% confluence because at each duplication the telomeres shorten, yielding to the onset of replicative senescence [61,62,63]. In the case of cancer cells, i.e., T, the telomeres shorten with doublings, and the cells will enter a state of crisis or death. Both normal and cancer cells could die whenever there is no space, because in vitro crowding situations lead some cells to undergo de-adherence and die [70]. Moreover, crowding-induced elimination facilitates the expansion of the malignant clones [78,79,80,81]. Due to the duration of the in vitro experiment (see Section 2.1) and for simplification, we discarded stem cells with ${\rho }_{max}=\infty$ in the model. The proliferation capacity parameter (${\rho }_{max}$) as attribute to T death shows that cells with a lower ${\rho }_{max}$ value tend to be at the edges of the clusters, probably competing for space with cells with a higher ${\rho }_{max}$ value (Figure 6e). As soon as the cells on the periphery die, the others proliferate and thus decrease their proliferation capacity. In this respect, the ignition of stem cells by space freed and/or their capacity to proliferate more than those with no stemness are in agreement with previously published theoretical models [64,65,66,67,68,69], experimental evidence [82], and the amalgamation of experimental data with mathematical modeling [83,84]. The proliferation capacity as a simpler rule to distinguish death from survival allowed us to observe a decrease in the H:T density ratio (from 10:1 to 4:1) (Figure 6g). It is biologically sound because pathways involved in contact inhibition might go through a transition towards the regulation of autophagy, which in turn can yield survival and proliferation [85]. Hitherto, we could only estimate ABM rules for proliferation and migration and several intrinsic conditions of cells beyond proliferation capacity play an important role in survival and/or death. Probably, a time series with more variables or data points than ours could display more basins of attraction; in this manner, the rules for survival and death could be estimated. In doing so, if one can estimate the rules by which cells proliferate, migrate, die, and survive, an implementation of life-like CA become possible [7]. In that scenario, it would be possible to map the ABM dynamic backward through algorithms like those presented at the software DDLAB (www.ddlab.org) in order to evaluate the dynamic upon or on the network [44,49] and to develop a bioinformatic pipeline to evaluate if an ABM is chaotic or ordered as well gain insight about the path to system attractors. Despite this, here we propose a possibility of describe TME regulations in terms of a logical rules from the tumor growth curve, moreover, the logical rules $(T\ge S)\vee (T\ge H)$ and $(T\le S)\vee (T\le H)$ serves to define how neighborhood state cells drive population dynamics. Either model-driven experimentation or data-driven modeling is a key way to improve and combine mathematical, biological, and clinical models [86]. Concomitantly, a discrete approach can also be combined with a continuous one [87,88], and it is possible to estimate important parameters that go beyond how dynamics change locally or globally. Our approach deserves room for further refinement because cell co-culture is an experimental system where many data points can be obtained, as well as images of spatial patterns. This makes it much easier to explore networks at the tissue scale, especially when it comes to figuring out how to reconstruct population dynamics.

5. Conclusion

In this work, we introduce a novel approach to cell proliferation dynamics. We present a approach to estimate ABM rules from the basins of attraction of a Boolean network and, also, align them with experimental conditions. By establishing a direct connection between our theoretical model and experimentally observed spatial patterns, we attain a significant milestone in the field of cellular analysis by effectively connecting computation and experimental methods. This modeling approach can make a significant advancement in our ability to comprehend and predict complex biological phenomena.