A computational model of the spatiotemporal adaptation of tumor cells metabolism in a growing spheroid

The Warburg effect, commonly depicted as an inherent metabolic trait of cancer in literature, is under intensive investigation to comprehend its origins. However, while the prolonged presence of excessive lactic acid production in tumors has been noted, it merely constitutes a fraction of the potential metabolic states cancer cells can adopt. This study aimed to elucidate the emergence of spatiotemporal diversity in tumor energy metabolism by expanding an existing model based on experimental facts. The resulting hybrid model integrates discrete formulations for individual cells and their processes, along with continuous elements for metabolism and the diffusion of crucial environmental substrates like oxygen, glucose, lactate, and the often underestimated acidity. This model enables simulation of a tumor spheroid, a standard experimental model, composed of numerous cells which can have distinct traits. By subjecting the spheroid to alterations of the environment such as cyclic hypoxia, acid shocks, or glucose deprivation, novel insights into metabolic regulation were obtained. The findings underscore the significance of the pyruvate-lactate interaction in governing tumor metabolic routes. Integrating acidity’s impact into the model, revealed its pivotal role in energy pathway regulation. Consequently, the conventional portrayal of a respiration/fermentation dichotomy proves inaccurate, as cells continuously and spatially adjust the ratio of these energy production modes, in contrast to abrupt, irreversible switches. Moreover, a cooperative cellular behavior akin to the reverse Warburg effect has emerged. This implies that the Warburg effect is not universally inherent to tumor metabolism, but a contextual, transient metabolic expression. Ultimately, the dynamic cellular-environment metabolic landscape influences cells’ survival under external conditions, with epigenetic regulations shaping their mobility potential within this landscape. While genetic mutations within tumor cells are undoubtedly present, this study shows they are not invariably essential for extreme metabolic modes or pathological characteristics to arise. Consequently, this research paves the way for innovative perspectives on metabolism, guiding tailored therapeutic strategies that consider not just patient-specific tissue attributes but also treat tumors as intricate ecosystems beyond their genetic diversity. Author Summary For years, scientists have been intrigued by the peculiar energy consumption patterns of cancer cells, such as the Warburg effect characterized by excessive lactic acid production. This study aimed to decipher the underlying reasons for the varying energy behaviors observed in different parts of tumors. Using a computational model, we simulated the collaborative dynamics of cells within tumors. The results revealed compelling insights. Two molecules, pyruvate and lactate, were identified as influential players in shaping energy utilization. Remarkably, the surrounding acidity was also found to exert a significant impact. Interestingly, tumor cells display a certain flexibility in their energy production strategies, adjusting according to prevailing conditions to maintain their survival and adaptability. Interestingly, cellular cooperation challenges the Warburg effect as an omnipresent phenomenon and reveals a transient nature. Our study underscores the significance of environmental influences, shedding light on the interplay between genetic modifications and the tumor environment in shaping cellular behavior. These findings hold promise for transforming cancer comprehension and devising treatments that tailor to both patients and the distinctive characteristics of their tumors.

1 Introduction Figure 1: Structure of the hybrid multiscale model. The simulation domain is defined by a grid in which the PDEs for the four environmental variables are solved. The tumor tissue is made up of interacting agents, the cells, with free positioning (lattice free). Each cell has its own characteristics and dynamics.
occurring during very frequent processes, on less frequent processes. For example, the phenotype of each cell 115 is evaluated only every 6 minutes of virtual time while their metabolism is evaluated every 1.5 seconds. In the 116 metabolism frame of reference, the phenotype therefore appears to be stable. Another very important point is the 117 saving in calculation time. If all the processes are evaluated at the smallest time scale, the computation time would 118 be too long without obtaining more realistic results. Hence, the following time spans were used for each modeled 119 biological process: 0.6 seconds for diffusion, 1.2 seconds for metabolism (resolution of differential equations = 2× 120 diffusion), 6 seconds for mechanics (10× diffusion) and 6 minutes for phenotype (600× diffusion). The substrates 123 In this study we focus on the cell energy metabolism. Therefore we chose to consider the three main following 124 substrates to fuel the cell energy: oxygen, glucose and lactate. We moreover include protons as the fourth en-  Non-homogeneous diffusion 155 In the Physicell version that we used (version 1.7.1), the solvers for diffusion equations were only considering a 156 constant diffusion coefficient. The initial simulations that we did, were not satisfactory as we could not generate The objective of this model was to look at the interaction between the gene regulatory network expressed with 177 ODEs and the metabolic expression. This model has more variables and is finer in the description of the reac-178 tions than our previous reduced model [4], but it keeps the same overall modeling philosophy. The parameters are 179 based on experimental data and the system as a whole is numerically stable. The system of differential equations  However, in its current version, the model was not fit to address some of our questions regarding spatiotemporal 192 metabolic heterogeneity. In particular, the time evolution of the different cellular configurations was not consid-193 ered. Instead, the probability of states via the topology of the attractors was privileged. Thus it was the derivative As a consequence, we took over the equations of this model by adding different aspects that were not covered 200 by the authors, specifically, the temporality, spatiality, adjustment of lactate transport mechanisms and influence of 201 the pH. in MATLAB code. It has been translated into symbolic code in Python, which is used to directly generate the C++ 205 code of the system of equations as well as the associated Jacobian matrix, for direct integration into PhysiCell.

207
The differential equations describe the one-or two-way reaction flows and take the following form [19] which here V max f and V maxr are respectively the maximum speeds of the reaction in the direct (forward) and reverse (back-211 ward) directions. K ms and K mp are the Michaelis constants of the substrate and the product respectively. and FADH 2 are each oxidized but do not produce the same amount of protons in the intermembrane space. The

248
NADH produces a little more which allows in fine to make the ATPase work and produce a little more ATP (ap-249 proximately 2.5 ATP for NADH and 1.5 ATP for FADH 2 ). In the model the FADH 2 is not taken into account, the 250 authors therefore establish, in a first reaction (r25), the account of the excess ATP produced by NADH (in relation 251 to FADH 2 , i.e., one molecule of ATP) thus creating the "complex II" mitochondrial unit. And, in a second step, they It is currently estimated that on average (again it depends on many conditions) for 12.5 ATP generated in the 257 mitochondria from a molecule of pyruvate, 80% comes from NADH (10 ATP) and 20% of FADH 2 (1.5 ATP). If we 258 reduce this calculation to the number of complex II formed, we therefore have: 1.5 AT P FADH 2 · 0.2 + 2.5 AT P NADH · 259 0.8 = 2.3 AT P/complex II in the electron transport chain. The final formula for the oxygen consumed can therefore be summarized as follows: That is to say :  • whatever the pH, the ratio of intracellular lactate/pyruvate concentrations is greater than approximately 10;

278
• the pyruvate concentration is kept relatively constant and remains above 0.1 mM;

279
• when the extracellular pH is acidic and the extracellular lactate concentration is higher than the intracellular 280 one or the lactate/pyruvate ratio falls too low, lactate is transported into the cell;

281
• when the extracellular pH is neutral and even if the extracellular lactate concentration is higher than the 282 intracellular one, very little lactate enters the cell.
288 Outward transport of lactate and proton r21 f orward = 1 e 2·(Lactate out −Lactate) + 1 · Lactate Lactate + 0.01 · Regulation lac/pyr with Lactate out the concentration of extracellular lactate and Lactate the intracellular one.

289
The first term of r21 f orward , makes it possible to reduce the outgoing flow when the extracellular lactate con-290 centration is close to or higher than the intracellular one and, conversely, to bring it to 1 when Lactate > Lactate out .

291
The second term cancels the first when the intracellular lactate is close to zero (so as not to artificially release lactate 292 which does not exist).

294
The last term manages the dynamic between lactate and pyruvate and results in: It is subtracted from 1 in order to increase the inflow when the concentration of extracellular lactate is close to or 317 greater than the intracellular one and conversely to bring it to 0 when Lactate > Lactate out . The third term cancels 318 the second when the extracellular lactate approaches 0 so as not to transport lactate inside when it is not present.

319
The last term Regulation pH e /pH i is defined by:

320
Regulation of the inward transport of lactate and protons by pH This last term varying between 0 and 1, promotes (the term approaches 1) the incoming transport of lactate and 322 protons when the difference between pH e and pH i weighted by pH eq is positive. Typically the value of pH eq is 1.1 323 which allows an equilibrium for an extracellular pH slightly more acidic than the intracellular one. black is slightly below its maximum speed of 1 when the extracellular lactate is close to 0 because the pyruvate is 326 very close to its minimum concentration and the lactate/pyruvate ratio is not large enough to compensate. Export 327 is therefore important due to a lack of lactate on the outside, but not maximum to avoid depleting all the pyruvate.

328
The pH being sufficiently acidic, the inflow (in red) tends towards its maximum speed of 0.2. For protons, the only fluxes having an impact are those of lactate transport (r21). As specified above, for an 334 excreted or integrated lactate molecule, a proton accompanies it. We would therefore tend to think that: In reality, this is not really the case because of the phenomenon of buffering (buffer effect of a solution). When 336 a small amount of proton is added to a solution, it is able to withstand a change in pH so close to an equilibrium Evolution of the concentration of extracellular H+ ions over time with Ka lac = 10 −3.86 mole/L, dissociation constant of lactate and Ka water = 10 −14 mole/L, dissociation constant of water. The constants of 1000 are here to convert molars into millimolars.

342
A note on the implementation in PhysiCell: when several cells are in the same place (above the same voxel), 343 they read the same concentrations at a time t and attempt to transport the same fluxes. In the event that there are 344 not enough nutrients for all the cells in the voxel, the risk is that the concentrations become locally negative, i.e.

345
the cells have ingested more nutrients than available at a moment t. To avoid this, a "voxel debt" system has been Evolution of the expression level of a gene over time The expression level of each gene X increases at a basal rate A of 0.005/min (normalized level/min) and is degraded at the same rate D. These values are calibrated on measurements carried out in previous studies and specified in the related publications. H(Y ) is a Hill function defining the regulation of the gene Y on the gene X: γ defines the type of regulation (positive > 1, negative < 1 and neutral =1), S the gene level Y at which the expression is half (similar to K m of a Michaelis Menten) and n the regulatory force (between 0.025 and 40).

357
All these equations are also defined in the Appendix S1. No changes have been made to these equations apart 358 for that of HIF. It is degraded at a constant rate in the original publication. We made it depend on the oxygen level 359 by adding the term H(O 2 ) to the HIF degradation term. is ATP. If the cell happens to be below a critical ATP threshold, it has a high probability of inducing necrotic death.

406
The more the level of ATP is below the threshold, the more this probability increases. This threshold is arbitrarily   Although the simulations can be made in 3D, we will only consider 2D simulations because of the high compu-439 tational cost. As a consequence, the cells can only move on the plane and the "2D spheroid" diameter will expand 440 faster than the real 3D one. The goal here is not to quantitatively compare the spheroids size and growth rate but 441 simply to find common emerging structures to appreciate their similarities. 449 Figure 6: Reference simulation of spheroid growth. Each image is spaced 50h from t = 0 to 400h. The color code is relative to each time step and indicates the age of the cells. The more purple a cell is, the older it is relative to other cells.

Gradients and radial profiles 450
The main advantage of a spheroid is that it is a simple tissue model with radial symmetry. Thus the majority of the 451 emergent heterogeneity is along this radial profile. This comes from the exposure to different gradients of substrates 452 that are differentially processed (produced or consumed) from the periphery to the core of the spheroid.
The first emerging profile is that of the distribution of cellular states (Fig.7). In this figure and on the spheroid 455 growth graph (Fig.8), the spheroid begins its growth with only proliferating cells. None entered the G0 phase and      The Warburg effect appears to dominate tumor cell phenotypes with few exceptions of the reverse Warburg 516 effect. This indicates that even cells previously labeled with respiration or pronounced respiration excrete lactate.

517
In fact, the latter secrete it in much smaller quantities than cells in hypoxia. In reality, it is not because lactate is

525
In Li and Wang [5] from which the model presented is derived, metabolic attraction zones are described but also 526 transition phases between these zones. Metabolic states are described as basins continuously linked to each other.

527
Four pools of attractions are identified (normal, oxidative cancer, glycolytic cancer, intermediate cancer) by the 528 expression levels of LDH and PDH genes. However, we here chose to replace normal by neutral so as not to force 529 a distinction between a healthy and cancerous mode. By choosing these two genes, it is the bifurcation at the level 530 of the fate of the pyruvate which is evaluated. The interest for looking at the levels of these genes in relation to the 531 reactions on which they act (respectively r12 and r18), is that the variation in their level of expression is "slower".

532
This makes it possible to look at the medium term effects of the environment on the metabolism of the cell, rather 533 than its immediate response to environmental fluctuations.

535
It is therefore interesting to look at the changes in the topology of this landscape, but mainly by looking at the capacity of the cells to move around, that is to say to reach the identified zones.  Cancer cells are often exposed to cycles of hypoxia, notably caused by a complex and constantly reconfiguring 554 vascular network [27]. Experiments to show the effect of these variations have been conducted with tissues in-555 termittently exposed to oxygen tensions ranging from 10 mmHg to 60 mmHg (some up to 160 mmHg but these 556 values are not physiological) [28]. These oscillations can have very fast (less than an hour), intermediate (tens of 557 hours) or slow (greater than a day) frequencies [28,27]. It is interesting to look at the repercussions of such oxygen 558 oscillations on metabolism and to look at the ability of cells to adapt more or less quickly to these changes. For   The phenotypes associated with these variations are now evaluated. Figure 15  Thus, fermentation is the dominant phenotype of the cell population in the hypoxic period with between 75% 599 (for the highest oscillation frequency) and 90% (for the lowest oscillation frequency) utilization. Conversely, the 600 Figure 15: Evolution of the proportion of metabolic phenotypes in the spheroid with or without hypoxia cycles. When a phenotype is described with a (-) it means that it is present but in a lower proportion than that marked with a (+).
"Fermentation (+) Respiration (-)" hybrid mode is predominant during the oxygenation period (between 40% and 601 55%). We can note a progressive fall in the use of respiration (majority at the start during the oxygenation period) 602 within the population whatever the simulation considered, due to the increase in the size of the tissue and the con- The 72-hour cycle in particular, creating longer fermentation times, generates more lactate but also has more  Here the pH has a double effect. By being identical everywhere, it initially allows cells to explore a wider 663 and more permissive metabolic zone. The reachability of these areas is then directly linked to access to the main 664 substrates, which are oxygen and glucose. Secondly, at a more acidic pH, it inhibits glycolysis. Whether this is 665 through some form of negative feedback (in the case where fermentation is responsible for the pH drop), or whether 666 it comes from an imposition of pH in the environment, the result is the same, fermentation is impaired. Figure 18: Evolution of the proportion of metabolic phenotypes in the spheroid population over time before and after the pH change in the whole environment.

Glucose depletion
In the same way as was done for pH, the following simulation aims to look at the consequences of glucose deple- The intrinsic heterogeneity of cell metabolism makes it difficult to derive a general law for tissue evolution in 745 response to environmental conditions. In particular, tissue dynamics quickly take precedence over the individual 746 metabolic dynamics of each cell. Thus, on many simulations, the behaviors that emerge are not always systemati-747 cally anticipated. However, from a large number of simulations that we realized the following trends emerge:  Table 2: General trends of simulations as a function of external conditions. The oxygen, glucose and lactate concentrations are reported qualitatively and the results column indicates the consequences of these conditions on the exposed cells.

749
Thus, lactate is generally excreted from the cell when oxygen is in low quantity in the environment or when 750 lactate is already present in the environment. Indeed, when lactate is present, the expression of the HIF-1α gene is 751 increased, which induces (although there is oxygen) an increase in glycolysis. If there is no lactate but glucose and 752 oxygen in sufficient concentrations, then the production is generally directed towards respiration, the production of 753 lactate is relatively low (but not non-existent).

755
On the other hand, if lactate is present in large quantities in the environment, the cells will tend to import it, 756 following the concentration gradient. However, this does not prevent the cell from continuing to produce it (by To illustrate this, a new simulation is presented. In this simulation, the oxygen is at atmospheric pressure with 768 160 mmHg i.e. in very great excess, whereas the glucose is only at 1 mM, which is a very small quantity. These (particularly in zone 1) is used to produce lactate and acidity which accumulate in these cells (zone 1 and 2) and in 791 the medium (at the balance between the two).
792 Figure 22: Cell phenotypes of the spheroid at 420h. On the left the cell states. In the center, the staining corresponds to cell phenotypes of fermentation and/or respiration. On the right the cells are stained in relation to their phenotype corresponding to the definitions associated with the Warburg effect. When a process is described with a (-) it means that it is present but in a lower proportion than that marked with a (+). The black cells are in necrosis.

793
In doing so they increase their level of pyruvate and more easily maintain their level of ATP. Surrounded by   : Flow of the reaction catalyzed by LDH (r12) and transport of lactate via MCT channels (r21) in the spheroid at 420h. The fluxes are expressed in mM/min. A negative value indicates that the reaction or transport is in the opposite direction (lactate to pyruvate and extracellular lactate to intracellular lactate). The few cells in green on the right spheroid represent those that produce much more lactate than the others (about 20x more). They appear in green so as not to interfere with the coloring of other cells.
These simulations show that it is uncertain to speak of the Warburg effect when thinking of targeting a single 811 type of metabolism in a particular way. If the Warburg effect appears on a macroscopic scale (at the tissue level) 812 this is not particularly the case at the level of individual cells. As we already explored this question in [1,2], the 813 definitions attached to it do not make it possible to specify the contours of a precise metabolic behavior at the 814 cellular level and even less of the threshold values (for example is the Warburg effect the production of X amount 815 of lactate in the presence of Y amount of oxygen). Even without defining these thresholds, taking the tissue as a 816 whole and simply considering lactate production (Otto Warburg's original observation) the Warburg effect is not 817 ubiquitous.

819
At the level of the metabolic landscape, three zones out of the four are identified (Fig.25). Initially (∼ 10h),  The simulations highlight the complexity of the multiple behaviors that can take place within the same tissue, 855 with cyclic hypoxia phenomena when the external conditions are at the limit barely sufficient to allow a few cells 856 to proliferate. There can also be phenomena of indirect cellular collaborations, of "metabolic symbiosis" [33], the 857 metabolic products of a cell exposed locally to certain environmental conditions benefiting another cell exposed to 858 different conditions.

860
As a consequence, it is possible to conclude that the metabolism is relatively robust in the face of environmental 861 disturbances, at least over short periods of time (<1 month). In reality, what introduces instability in the metabolic 862 organization of a tissue does not come intrinsically from the metabolism, but from the existence of abnormal con-863 ditions in which the cells are pushed (the depletion of glucose in the environment resulting in death is an example).

864
A spheroid, too, is an abnormal tissue configuration due to its high density or the absence of vascularization. The