Physics-based tissue simulator to model multicellular systems: A study of liver regeneration and hepatocellular carcinoma recurrence

We present a multiagent-based model that captures the interactions between different types of cells with their microenvironment, and enables the analysis of the emergent global behavior during tissue regeneration and tumor development. Using this model, we are able to reproduce the temporal dynamics of regular healthy cells and cancer cells, as well as the evolution of their three-dimensional spatial distributions. By tuning the system with the characteristics of the individual patients, our model reproduces a variety of spatial patterns of tissue regeneration and tumor growth, resembling those found in clinical imaging or biopsies. In order to calibrate and validate our model we study the process of liver regeneration after surgical hepatectomy in different degrees. In the clinical context, our model is able to predict the recurrence of a hepatocellular carcinoma after a 70% partial hepatectomy. The outcomes of our simulations are in agreement with experimental and clinical observations. By fitting the model parameters to specific patient factors, it might well become a useful platform for hypotheses testing in treatments protocols.

in different degrees. In the clinical context, our model is able to predict the recurrence of a hepatocellular carcinoma after a 70% partial hepatectomy. The outcomes of our simulations are in agreement with experimental and clinical observations. By fitting the model parameters to specific patient factors, it might well become a useful platform for hypotheses testing in treatments protocols.

Author summary
We introduce an off-lattice agent-based model to simulate tissue-scale features that emerge from basic biological and biophysical cell processes. In order to calibrate and validate our model, we have considered the liver regeneration response after a 30% partial hepatectomy in which the liver recovers its original volume due to the hypertrophy of the hepatocytes. Subsequently, we have modeled the same process but after a 70% partial hepatectomy, in which the liver recovers its original volume due to the hypertrophy and the proliferation of the hepatocytes. Unfortunately, the precise mechanisms of initiating, promoting and terminating regenerative responses remain unknown. As a consequence, we have proposed a modeling approach in which such processes are regulated by a hypothetical substrate that diffuses in the cell microenvironment. As a further test, we have, in one hand, implemented our model to predict the liver response after a 50% partial hepatectomy and, on the other hand, explored our model's ability to account for the recurrence of a hepatocellular carcinoma. The outcomes of our simulations agree with experimental data and clinical observations, which comes to underline the significant descriptive and predictive power of this computational approach. Even though our model needs to be further extended to incorporate patient specific clinical data, these results are a promising step in the direction of a personalized estimation of tissue dynamics from a limited number of measurements carried out at diagnosis.

Introduction
Many significant multicellular system problems such as tissue engineering, evolution in bacterial colonies, and tumor metastasis can only be understood by studying how individual cells grow, divide, die, move, and interact. Tissue-scale dynamics emerges as cells are influenced by biochemical and biophysical signals in their microenvironment, bearing in mind that cells themselves continually remodel their own microenvironment. Thus, the ideal scenario to study a multicellular system's biology must simultaneously address: tissue microenvironments with multiple diffusing chemical signals (e.g., oxygen, drugs, and signaling factors), and the dynamics of many mechanically and biochemically interacting cells.
To that aim, mechanistic dynamical systems models described by ordinary differential equations have been developed [1][2][3][4][5][6][7][8][9][10] (cf. those reviewed in [11] and [12]). Even though these models are very useful, they lack the spatial resolution that would enable the examination of intratumoral heterogeneity and its correlation with treatment efficacy. This is a relevant feature since intratumoral heterogeneity has become a central element for understanding important aspects of cancer treatment such as drug resistance and biomarkers [13,14].
A widely used modeling paradigm in the study of complex biological systems is the agent-based model (ABM) [15,16]. ABMs are implemented mainly to simulate the actions, behaviors and interactions of autonomous individual or collective entities, with the aim of exploring the impact of an agent or a type of behavior in the system. Agents are independent units trying to accomplish a set of goals. The purpose of an ABM is to explore variations in the system behavior due to agent characteristics or rules. These attempt to emulate the general behavior of the system and predict the patterns of complex phenomena. Agents behave independently, but they react to the environment, modify system properties, and incorporate new agents. They also have the ability to "learn", that is, to avoid previously unsuccessful decisions and favor successful ones, as well as to "adapt", i.e. change their behavior in response to variations of the properties of the system. Their basic advantage is that they provide predictive power on a large scale [17].
In many cases the system being modelled is usually comprised of millions of components. In those cases a great level of abstraction/simplification must be applied to the model to render it useful [18] (we refer the reader to [19] for as extensive discussion on the accuracy of computational models).
For that reason, once a preliminary model has been constructed it must be subjected to verification, calibration and validation. Verification is the process of determining how accurately prior knowledge and underlying assumptions have been translated into mathematical form. Calibration is the process by which parameters in a model are adjusted so as to match model performance to experimental data. Finally, model validation is the process of evaluating model performance against the primary design goal. In the case of biological models, this usually aims at achieving a close match between model and experiment [20].
In a medical context, ABMs should allow simulating clinical trials in sufficient detail to study the subject's response to changes in therapy in simulations, rather than in patients. ABM have been used to study many different aspects of cancer, including tumor growth, cancer cell migration, metabolism, evolutionary dynamics, metastasis, angiogenesis and the role of cancer stem cells [21][22][23][24][25][26][27][28][29][30][31][32][33][34]. The reader is referred the extensive review by Norton and coworkers [35] for a more detailed presentation of the potential of ABMs in the context of cancer modeling.
Here, we will present a mechanistic off-lattice agent-based model to simulate tissue-scale features that emerge from basic biological and biophysical cell processes. This ABM will be validated and put to test for modeling the unusual ability of the liver to regenerate [36]. Even when 70% of its mass is surgically removed, the remnant portion expands to compensate for the lost tissue and functions [37]. The multilobular structure of the liver allows the surgical resection of a lobe of choice to achieve different degrees by partial hepatectomy (PH). Because the resection of lobes does not induce damage to the remaining tissue, PH is a clean model. Therefore, liver regeneration after PH has long been an excellent experimental model for tissue regeneration. Furthermore, although the liver consists of various types of cells, hepatocytes account for about 80% of liver weight and about 70% of all liver cells [38]. Thus, hepatocytes provide an ideal starting point to study the relation of organ size with number and size of cells.
It has been generally accepted that liver regeneration depends mainly on the proliferation of hepatocytes [37,39,40]. However, there are several reports showing hypertrophy of hepatocytes in the regenerated liver [41][42][43]. Miyaoka et al. [44,45] performed a series of experiments and found that although a number of studies indicated that almost all hepatocytes proliferate after 70% PH, cellular hypertrophy significantly contributes to liver regeneration as well. They showed that hepatocytes undergo cell division only about 0.7 times on average in the regeneration from 70% PH, and that at early stages, the regeneration totally depends on the hypertrophy of hepatocytes.
In contrast, liver regeneration after 30% PH is solely due to hypertrophy. Therefore, liver regeneration process is a perfect scenario to test and calibrate our model. On the other hand, post-hepatectomy liver failure is a serious complication after liver resection and its incidence varies from 1.2 to 32% [46][47][48][49]. It is defined as functional deterioration of the liver associated with an increased international normalized ratio (INR) -a measurement of how long it takes blood to form a clot-, and hyperbilirubinemia typically starting after the fifth postoperative day [46]. There are recommendations that post-hepatectomy liver failure could be prevented if the remnant liver size is above 20% of its original size in patients with normal liver function, and 30 − 40% in patients with steatohepatitis or cirrhosis [50,51]. Nevertheless, even with adequate preoperative assessments, post-hepatectomy liver failure is a major contributor to mortality rates of up to 5% after liver resection [52,53]. Various patient-related factors (age, comorbidities such previous chemotherapy, cirrhosis, fibrosis, cholestasis, and steatosis), and surgery-related factors (extent of resection, blood loss, and ischemia reperfusion injury) affect the regenerative capacity of the remnant liver [54,55]. Given all these numerous factors, estimating the adequate extent of the hepatectomy, and individual regenerative capacity, remain significant challenges for clinicians and scientists. It is to be stressed here that the regeneration process is controlled by different transcriptomic signatures depending on the intensity, duration and etiology of liver injury [56].
Consequently, different transcripts can modulate the results of the regeneration processes.
In any case, according to the Barcelona Clinic Liver Cancer (BCLC) staging system [57], hepatic resection (partial hepatectomy) can be considered as a curative treatment for patients with stage 0-A uninodular hepatocellular carcinoma (HCC) who maintain preserved liver function and without portal hypertension. The prognosis of HCC patients has improved because of advances in radiologic assessment, patient selection, operative techniques, and perioperative care [58,59]. On the other hand, long-term prognosis of patients with HCC after liver resection is often affected by high tumor recurrence rates that reach 40 − 70% within 5 years [60]. This is an issue that must urgently be addressed, and where we believe that a well calibrated computational model that succeeds on modeling the liver regeneration processes and HCC evolution, would complement the clinical trials, once fed-in with specific patient data. Moreover, liver regeneration is the basic element for the maintenance of liver function and size during homeostasis and liver injury (acute and chronic).
Understanding the mechanisms of hepatic regeneration, from its cellular origins and signaling mechanisms, is essential to design specific regeneration models. Making these models available in a more "accessible" way that allows for a quicker evaluation of the influence of certain factors on regeneration processes, unlike the "classical" models, would then facilitate the optimization the different diagnostic strategies (time intervals for hepatocarcinoma screening) and treatments (recurrence rate compared to ablative therapies, immunotherapy, response to sequential therapies). The approach here presented might well constitute a potential tool to evaluate biomarkers, such as circulating

Biological model
The liver is a highly complex organ, which removes drugs and toxins from the blood. It is characterized by its multi-scale architecture (figure 1) which consists of four lobes: the right lobe, the left lobe, the caudate lobe, and the quadrate lobe, which are further divided into eight segments based in the Couinaud system, also known as hepatic segments [65]. Each segment has its own vascular inflow, outflow and biliary drainage. The division of the liver into independent units means that segments can be resected without damaging the remaining segments. Hepatic parenchyma is organized in repetitive functional units called liver lobules. The lobules are roughly hexagonal prisms, and consist of plates of hepatocytes, and sinusoids radiating from a central vein towards an imaginary perimeter of interlobular portal triads. The portal triad is a distinctive component of a lobule, which can be found running along each of the lobule's corners. It consists of the hepatic artery, the portal vein and the common bile duct. The central vein joins the hepatic artery and the portal vein to carry blood out from the liver [66]. A recent work reported that this particular structures remain during the liver regeneration process [67]. Even though hepatocytes (liver parenchymal cells) account for about 80% of liver weight and about 70% of all liver cells, the liver also has other type of cells named nonparenchymal cells: endothelial cells, Kupffer cells (macrophages resident in the liver), and biliary-duct cells [38].
One of the main characteristics of the liver is its high regenerative capacity after injury. Even when 70% of its mass is surgically removed, the remnant portion expands to compensate for the lost tissue and functions [37,68]. Liver resection is the most common liver surgery and consists of removal of liver tissue due to focal lesions, most often malignant tumors and living liver donation [69]. The multilobular structure of the liver not only allows the surgical resection of a lobe of choice to achieve different degrees by partial hepatectomy but also the resection of lobes does not induce damage to the remaining tissue. The extent of resection is determined by the size and location of the focal lesion and the estimated function of the future liver remnant [70]. Prior to liver resection, surgeons have to assess the patient's individual risk for postoperative liver dysfunction. In case of malignant tumors, surgeons have to identify the surgical strategy best suited to allow radical oncological removal in order to avoid recurrence, without putting the patient at risk of postoperative liver failure due to excessive removal of liver mass [71][72][73].
The liver regenerates in a highly organized fashion after surgery [68]. The human body responds to partial hepatectomy not by regenerating lost segments but by inducing hyperplasia in the remnant liver [36,74,75]. The anatomical structures of a liver that has undergone partial hepatectomy are therefore distinctly different from those of the original liver. The process of restoration of liver volume is initiated by the replication of various types of intrahepatic cells, followed by an increase in cell size. Nonparenchymal cells replicate in a delayed fashion. After replication has Figure 1: Schematic representation of the multi-scale liver architecture. The human liver is divided grossly into four parts or lobes. The four lobes are the right lobe, the left lobe, the caudate lobe, and the quadrate lobe. Seen from the front the liver is divided into two lobes: the right lobe and the left lobe. It is further divided in eight functionally independent segments based in the the Couinaud classification of liver anatomy. At the microscopic (histological) scale, the liver is organized in repetitive functional units called liver lobules, which take the shape of polygonal prisms (typically hexagonal in cross section). Each lobule is mainly constituted by hepatocytes and it is centered on a branch of the hepatic vein called the central vein which is interconnected with the interlobular portal triads: the hepatic artery (red), the portal vein (blue), bile duct (green). been completed, growth consisting of an increase in cell size occurs over several additional days.
The initiation and synchronization of replication in different types of hepatic cells depend on the extent of the resection, tissue damage, or both. Low-grade tissue damage (e.g., toxic or ischemic injury) or a relatively small resection (removal of less than 30% of the liver) substantially reduces the replication rate, which also appears to be less synchronized than after a large resection (removal of 70% of the liver) [36,68,75].
Since all our results are compared with experimental data obtained from rodent models, it is important to mention that the process of hepatic regeneration in rodents and humans is similar, and the results obtained from rodents could be applicable to the human liver [76]. Moreover, the rat liver architecture can be compared with the human liver segmentation defined by Couinaud [77].

Computational model
Our model is implemented resorting to an object oriented programming model, and to that aim we have used C++11 language. Simulation CPU time depends on model parameters such as domain (lattice) size, cell number and simulation length (in time); a typical simulation run takes approximately 6 h on a single core of an Intel i7-10510U CPU. Model visualization is performed with Ovito [78], Paraview [79] and Matplotlib [80].
In order to reduce the computational burden of the simulations, an abstraction process was necessary to go from the biological to the computational model. First, we disregard the explicit liver geometry, instead we have chosen a reduced spherical model. This simplification is possible because, as it was mentioned in section 3.1, after a PH the liver does not regenerate the lost segments, i.e. does not recover the original shape. Instead it just recovers its original volume by hyperplasia of the remaining lobes. Subsequently, we rather focus our attention on the liver parenchyma instead of liver lobes. In our computational model, hepatic lobules are hexagonal prisms delimited by an imaginary perimeter of interlobular portal triads. The central vein that carries the blood out from the liver as well as the liver sinusoids are not explicitly modeled as the portal triads are. We rather model this behavior by tuning the effective oxygen difussion and decay constants. A discussion on the implications of these simplifications on the results will be presented in section 5.
In the following subsections we will describe the methods implemented to model diffusion and cellular mechanics, as well as the mathematical models to predict tissue growth kinetics. For further details, we refer the reader to the supplementary material. A schematic representation of the inner workings of our model is depicted in Fig. 2. The key elements of the model are described in what follows. Figure 2: Main loop flow diagram. Blue box represents the start of the program. Red box represent the diffusion processes. Green box and orange box describe the cell mechanics and cycling processes respectively. Finally, yellow boxes represent the data saving process. After initializing the microenvironment, the cells, and the current simulation time t = 0, our model tracks (internally) t mech (the next time at which cell mechanics methods are run), t cycle (the next time at which cell processes are run), and t save (the simulation data output time), with output frequency ∆t save . % represents the modulo operation.

Diffusion solver
We model the diffusion of chemical substrates in the tumor microenvironment as a vector of reaction-diffusion partial differential equations for a vector of chemical substrates, ρ. It is discretized over a Cartesian mesh for computational convenience, in such a way that each voxel (volumetric pixel) stores a vector of chemical substrates. Each substrate diffuses and decays, and can be secreted or uptaken by individual cells at their specific positions.
We use a first order, implicit (and stable) operator splitting, allowing us to create separate, optimized solvers for the diffusion-decay, cell-based source/sinks [81]. The diffusion-decay terms are solved using the finite volume method [82], further accelerated by an additional first-order spatial splitting in the x−, y− and z−directions via the locally one-dimensional method (LOD) [81,83]. For each dimension, we solve the resulting tridiagonal linear systems with the Thomas algorithm [84].
We also implement the so-called Dirichlet nodes, so that substrate values at any voxel within the simulation domain can be overwritten to turn the voxel into a continuous source of substrates. This is particularly useful to model the effect of blood vessels, or when applying Dirichlet boundary conditions.
For computational efficiency we use thread parallelization to relevant loops, e.g. in many instances of the Thomas solver when solving x−diffusion across multiple strips of the solution domain.
This methods were already implemented and successfully tested by Ghaffarizade et al. [85], therefore, we have validated the numerical accuracy of our solver by comparing our results with those found in Reference [85]. For further details, please refer to the supplementary material (S1 Text)

Cell agents
Since we are implementing an agent-based model programmed in the context of an object oriented approach, each cell is an agent implemented as a software object that acts independently.
Like most classes in software it has attributes, i.e. its own internal variables that each specific agent is allowed to manipulate on its own. It also has methods, which are functions that act upon the attributes. In the context of cell biology, relevant attributes might be: position, size, cell state regarding the cell cycle, etc. The cell cycle is an object with the aforementioned attributes. The cell class have methods that represent cellular processes such as, death, growth and volume change, and are coordinated by the cell cycle object.
One of the main features of our model is that cells are off-lattice. Consequently, they are not confined to a particular lattice or spatial arrangement, they move in a continuous fashion through all space positions, and therefore underlying possible artifacts associated with the chosen lattice structure and spacing are removed.
Based on previous cell-based models [33,86,87], we have modeled the cell behavior as follows: Cell cycle: We model the cell cycle as a directed graph, in which the nodes represent the phases and the edges the transition rates between them. These transition rates can take stochastic or constant values. Moreover, any of the cell cycle time scales can be adjusted at the beginning of the simulation to match different types of growth and they can also be adjusted at any time on an individual cell in order to reflect its microenvironment influences.
Our model allows to implement different types of cell cycles based on different parameters.
Following Miyaoka and coworkers [44,45], we can base our cell cycle on a tracking of the expression of protein Ki-67 [88], involved in cell proliferation. This tracking is performed thanks to a nuclear stain of Ki-67.
Finally, since at cell scale death is not an instantaneous event but a process, we model death cycles such as Necrosis and Apoptosis as we did with cell cycles, by using directed graphs, death is part of the cell cycle after all. The entry to the death cycles depends on the resources, for example oxygen, drugs, therapies, etc.
Cell volume: To model cell volume variation, each cell tracks V (total volume), V F (total fluid volume), V S (total solid volume), V N S (nuclear solid volume), V CS (cytoplasmic solid volume), V N (total nuclear volume), and V C (total cytoplasmic volume). Key parameters include nuclear solid, cytoplasmic solid, and fluid rate change parameters (r N , r C , and r F ), the cell's "target" fluid fraction f F , target solid volume V * N S , and target cytoplasmic to nuclear volume ratio f CN . For each cell, these volumes are modeled with a system of ordinary differential equations that allow cells to grow or shrink towards a target volume. These parameters are updated as the cell progresses through its current cycle or death cycle.
Cell mechanics: In our model, as in [33], we use the piecewise polynomial force which is constructed as the sum of a positive adhesive and a negative repulsive polynomial force contributions.
It is important to note that repulsive forces are really an elastic resistance to deformation. To compute these forces we use adhesive and repulsive interaction potentials functions.
We assume that three types of forces act upon each cell. First we have a drag force, which represents dissipative drag-like forces such as fluid drag and cell-extra cellular matrix adhesion. We then have neighboring cell mechanical forces. In the simplest case these involve repulsive forces due to limited cell compressibility, but they usually also include cell-cell adhesion. We use interaction potentials that depend upon each cell's size, maximum adhesion distance, adhesion and repulsion parameters and distance to other cells. Finally, the third force acting on the cells is the cell-basement membrane interaction.
Since the cell microenvironment has a very low Reynolds number (the Reynolds number describes the ratio of inertial to viscous forces) [89], inertial effects such as acceleration are neglected. This is commonly known as the inertialess assumption, and implies that forces equilibrate at relatively fast time scales in contrast to the processes involved in cell cycling, death cycling, and volume variation.
We finally use the Adams-Bashforth method [90] to solve the mechanics equation to enhance computational efficiency.
We refer the reader to the supplementary material (S1 Text) for further details on the implementation of the piecewise polynomial force model.
Cell secretion and uptake: This is one of the most important data structures of the cell because it links the cell with its microenvironment. We solve a vector of partial differential equations, which reduce to the addition of a cellular secretion/uptake term to the diffusion equation described in This is very important since most of the cellular processes depend on the substrates that diffuse in the microenvironment. For example, it is well accepted that after a partial hepatectomy, the liver undergoes to cytokine-and growth factor-mediated regeneration processes [91]. However, most of the mechanisms of initiating and promoting regenerative responses as well as the termination of liver regeneration remain unknown [44,45]. In this work cellular proliferation is assumed to be controlled by a growth factor that diffuses through the microenvironment. This growth factor is only considered as an abstract parameter which encompasses all the underlying molecular mechanisms involved in the liver regeneration process. The cell cycle entry rate is proportional to this factor in the following way: where t K − is the mean time each cell spends in the non-proliferative phase (see section 3.2.2), which can be experimentally monitored using the Ki-67 cellular marker. GF is the current growth factor concentration value in the cell's voxel, GF prol is the proliferation threshold, i.e. the growth factor value below which the proliferation ceases and GF * is the proliferation saturation value, above which the proliferation rate is maximized. Therefore, based on the growth factor concentration, the hepatocyte will enter either the hypertrophy phase or the proliferation phase. A similar approach can be done based on the oxygen concentration, however, instead of influence on the decision about whether or not to proliferate, oxygen will accelerate the phase entry. Please refer to the supplementary material (S1 Text) for further information about how cell cycles and transition rates will be modified based on chemical substrates concentrations.

Growth estimations
To complement our model in the prediction of tumor growth and/or tissue regeneration we use a mathematical model known as the Gompertz model [92][93][94][95]. This model assumes initial exponential growth, but as the tumor grows the volume-doubling time increases due to lack of nutrients and subsequent cell death, by which the growth rate shifts towards linear regime, finally reaching a plateau [92]. This is given by: where the parameter K is the carrying capacity of the tissue, which is the highest possible volume, V 0 is the initial volume of the tissue and α is the specific growth rate [96,97] which is defined by Here V 1 and V 2 are tumor volumes at the measure times t 1 and t 2 respectively. This parameter α determines how fast the tumor reaches the carrying capacity K and its measured in (%/days).

Results
We first attempted to define a baseline scenario of the liver physiology. As it was mentioned before, we focus our attention in the liver parenchyma which, in our model, is made up of hepatic lobules which are hexagonal prisms delimited by an imaginary perimeter of interlobular portal triads. We idealized this dynamic vasculature architecture (figure 3a) by using the Dirichlet nodes as shown in figures 3b and 3c, and defining the distance between triads on micrograph analysis [66]. Figure 3b shows a transversal cut of our simulated liver in which we can observe the hepatic lobules architecture. Blue dotted lines were drawn just as a guide to the eye. Pink spheres represent the hepatocytes and the white squares represent the portal triads that oxygenate the tissue. Figure 3c shows the heat map of the oxygen diffusion in the liver micoenvironment. The oxygen diffuses from the portal triads (Dirichlet nodes), there is no diffusion from the boundaries of the simulation box. The grid we used to simulate liver regeneration is 1 mm 3 in total size, orders of magnitude smaller than an actual liver. While the model does not place limitations on the liver size, the implementation is obviously constrained by the size of the computer memory. However, this limitation can be to some extent circumvented by considering the sample region as representative of a significant subregion of the liver. Obviously size effects must be investigated in further detail to the extent that new algorithmic and hardware improvements become available (e.g. by exploiting massively parallel CUDA [98] cores and/or Tensor cores in new GPU/TPU architectures).
All the parameters values that were used in the following examples, and the corresponding sources are listed in section 4 of the supplementary material (S1 Text).

Liver regeneration
Our first experiments aim at assessing the ability of our model to describe the dynamics of the liver regeneration process. We performed the in silico version of the experiments of Miyaoka et al. [44,45], with the novelty that hepatocytes secrete and are sensitive to a growth factor that diffuses through the microenvironment creating a gradient. We assume that when the liver is intact, the microenvironment is in homeostasis and, although they are metabolically active, hepatocytes remain dormant in the cell cycle. When a partial hepatectomy occurs, each cell secretes a constant amount of growth factor depending on the extent of the injury. They will undergo either hypertrophy or proliferation as regulated by the concentration of the growth factor. This is illustrated in figure   4. In panel (a), we can see a peak of growth factor when the partial hepatectomy occurs. A 30% PH does not generate enough growth factor to make hepatocytes proliferate (figure 4b). On the contrary, a 70% PH not only produces enough growth factor to make the hepatocytes proliferate, they also undergo an hypertrophy process (figure 4c). Finally, as the liver regenerates, this growth factor decreases.

30% partial hepatectomy
We first analyzed liver regeneration after a 30% PH. We have set the hepatocytes cell cycle time 33.6 hr as stated in ref. [99] but as mentioned in section 3.2.2, this could vary depending on the oxygenation of the tissue and the growth factor concentration. Figures 5a and 5b illustrate the regeneration process in a qualitative and a quantitative manner respectively. We found that liver volume increased from 1 day after the PH and reached a plateau of 0.93-fold of the liver initial An animation of the liver regeneration process after a 30% partial hepatectomy can be seen in the Supplementary Material S1 Video.

70% partial hepatectomy
After a 70% PH, we found that liver volume increased from day 1, reaching a plateau with a total increase of 0.72-fold of the liver initial volume in day 7. Figures 6a and 6b illustrate this process in a qualitative and a quantitative fashion respectively. Hepatocytes entered the cell cycle 2 days after the PH, as shown in figure 6e. Although there were a few active proliferating hepatocytes on day 1, liver volume had increased considerably by that time. Then we measured the area of the hepatocytes (figure 6c) and found that hepatocytes increased their volume by 2.0-fold the first Figure 5: 30% PH (a) Qualitative representation of the liver regeneration process after a 30% PH.
(b) Fold-increase in the liver size. Observational data is shown in gray and represents the weight of the liver measured by Miyaoka et. al [44,45]. Simulated data is shown in red and represents the liver volume. (c) Quantification of the hepatocytes area during liver regeneration. day and gradually decrease by 1.5-fold the 14 th day after the PH (figure 6d). These results are in agreement with those obtained by Miyaoka et al. [44,45], and indicate that proliferation of hepatocytes alone could not account for the recovery of liver mass after a 70% PH.
The animation of the liver regeneration process after a 70% partial hepatectomy, can be seen in the S2 Video.
Although the exact mechanisms underlying liver regeneration have not yet been fully characterized, studies have shown that after 70% PH many of the upregulated growth factors in a regenerating liver are known for their angiogenic properties in vivo. For instance, vascular endothelial growth factor (VEGF) is upregulated after PH [100][101][102]. It is a major pro-angiogenic factor [103] and is thought to improve sinusoid reconstruction during the liver regeneration process [104]. In the processes of angiogenesis associated with tissue regeneration, two phases are described; an induction phase and another proliferative angiogenesis phase. In post-hepatectomy liver regeneration models, the first is calculated during the first 4 days and the second between days 4 − 12 [105]. We have quantified blood vessel formation by counting the Dirichlet nodes added during the liver regeneration process. Figures 7a and 7b represent this process in a qualitative and a quantitative way, respectively. On one hand we found that during the regeneration process, blood vessels keep the hepatic lobules architecture (figure 7a). To model this behavior, we have labeled each voxel according to its potential to become a blood vessel, inspired by the hepatic lobule architecture.
During the simulation, if the voxel is tagged as a potential blood vessel and contains five or more cells, it turns into a Dirichlet node, that will provide oxygen to the cell microenvironment. As shown in figure 7b, we found that the number of Dirichlet nodes increased significantly during the first 3 days following PH until they finally reach a plateau. These results are in agreement with ref. [106], in which it was found that the micro-vasculature density increased significantly during the first 3 days in mice subjected to PH. This proves that our model has the ability to simulate cell population behavior, and subsequently it also can predict blood vessel formation.

Test of a 50% partial hepatectomy
Although 70% PH is the most studied instance of liver regeneration [36,107], resection of approximately half volume of the donor liver is more common in the case of living donor liver transplantation setting. [68]. Moreover, as it was mentioned before, resection must be more conservative in the presence of underlying liver diseases or in elderly patients (e.g., ≥ 70 years of age) [68]. Major (> 50%) hepatectomy in the presence of cirrhosis or steatosis significantly increased morbidity [108]. With the presence of liver steatosis, 30% or more of the remnant liver should remain in order to maintain viability. Furthermore, studies have revealed an increased benefit of ischemic preconditioning in patients with hepatic steatosis who had lower resected liver volume (< 50%) [109], and extensive resections are generally not recommended for patients with cirrhosis [110]. Hemihepatectomy (i.e. 50% PH) has now been successfully and frequently used for surgical removal of liver associate tumors and cancers [111,112]. However, to the best of our knowledge, a study of liver regeneration dynamics after a 50% PH, such as the one presented in [44,45] for 30% and 70% PH, has not been reported yet. Using the calibration of the previously studied PH, we have applied our model to study the liver regeneration process after a 50% PH.
We found that the liver volume increased from day 1 until it reaches a plateau with a total increase of 0.86-fold of the liver initial volume, within 5.5 days. Yoshioka et al. [113] reported that 3 days after a 50% PH, the remnant liver weight reached 0.72-fold ±0.05-fold of its original calculated weight, which is similar to the 0.78-fold that our model predicted for the third day after the PH. As shown in figure 8a, we can apply a polynomial fitting that matches our simulations and estimates the outcomes of the hepatectomies. Figure 8b shows that, similar to the 70% PH, the hepatocytes not only undergo hypertrophy, they also proliferate. We found that hepatocytes  It is important to mention that the model is calibrated based on an normal hepatic parenchyma, i.e. without liver cirrhosis, with preserved liver function and no signs of portal hypertension, therefore its predictions correspond to that specific case. However it could be adjusted to model the presence of underlying liver diseases.

Recurrence of hepatocellular carcinoma
With the scenario of the liver regeneration working properly, we are now interested in study tumor recurrence. The need for extended liver resection is increasing due to the growing incidence of liver tumors in aging societies [71], however, the resected volume not only depends on the tumor volume itself but also on the patient's liver overall health. Prior the resection, surgeons not only have to assess tumor resection to avoid recurrence but also the patient's individual risk for postoperative liver dysfunction. It is well known that planning for a safe resection of a liver tumor with a large future liver remnant reduces the risk for postoperative liver failure but increases the risk of tumor recurrence. In contrast, planning for an oncologic radical surgery requires to remove a safety margin. Extending the safety margin in case of a centrally located tumor leads to a substantially extended resection leaving a rather small future liver remnant behind, which increases the risk of postoperative liver failure [71]. In a patient with uninodular hepatocarcinoma or up to 3 nodules smaller than 3 cm each without macrovascular infiltration, without distant metastasis (BCLC Stage 0-A), with preserved liver function and without portal hypertension, a liver resection of up to 60 − 70% is a feasible scenario considered in the BCLC classification [57], as long as it presents a normal hepatic parenchyma. Preexisting liver disease such as steatosis increases the risk for postoperative liver failure and might therefore call for a smaller PH compared to livers without preexisting diseases, however, the study of liver and hepatocellular carcinoma dynamics with preexisting diseases exceeds the scope of this paper.
Following the previous case in which we have considered a normal hepatic parenchyma, we have seeded a remaining tumor clone in order to model a recurrence of a hepatocellular carcinoma after a extended resection, i.e. a 70% PH. As proof of concept, we have assumed that preserved liver function without portal hypertension could approximate the behavior of a liver with preserved or normal parenchyma.
Cancer cell cycle has a duration of 38.6 hr [99] but the rate of cycle entry scales proportionally to oxygen concentration. As a first step, we have randomly seeded the tumor clone over the liver surface in order to test if its initial location would change the simulation outcomes. We have performed 40 simulations and obtained a 59% standard deviation of the tumor final size. Figure 9 shows the smallest (a) and the biggest (b) tumors, 30 days after the PH. Hepatocytes are drawn transparent to have a better view of the blood vessels (red tubes) and the tumor (blue cells) growth. We observe that the tumor volume vary depending on the location of the residual tumor clone, whether it is in the periphery or the center of the liver surface. This might be due to the process of blood vessels generation implemented in our model (sec. 4.1.2). When the tumor clone is located in the center of the liver surrounded by hepatocytes, it will have to compete for oxygen with the surrounding cells (i.e. share the preexisting blood vessels), reducing the cell cycle entry. On the contrary, when the tumor clone is located in the periphery, it will grow outwards. That means that the tumor cell not only has more oxygen for itself due to the liver blood vessels, but as the tumor gets bigger, it generates its own blood vessels, as shown in figure 9b, which increases the cell cycle entry. In clinical practice, among the factors that have been reported as being associated with early or late recurrence, there are characteristics such as tumor size and number of tumors before PH, micro and macrovascular invasion, degree of tumor differentiation, alpha-fetoprotein levels [114][115][116][117][118][119][120] but how tumor location in the liver will influence the HCC recurrence has been poorly described and with heterogeneous results [121][122][123][124].
We have considered the mean value of our simulations and analyze the recurrence dynamics.
We observe that, despite the initial location, there is a delay in the cancer cell reactivation in comparison to the liver regeneration process. Figure 10a shows the process in a qualitative fashion.
Cancer cells start growing after the liver finishes its regeneration process. In this case, they grow inwards, following the increase of oxygen concentration, and consequently they proliferate towards the blood vessels (the animation of the HCC recurrence can be seen at the S3 Video.). The growth of the tumor cells showed in figure 10b, allowed us to compute the specific growth rate of the tumor (eq. 3), which give us α = 0.053 %/day. By using the Gompertz Model (eq. 2), we can predict the  [126][127][128], in which the earliest time from surgery to first recurrence was 45 days. It is considered an early recurrence.
Our model has thus proven its ability to estimate the growth kinetics of the tumor based on its early stage growth. It could well turn into a useful tool to determine the optimal follow-up interval after the PH. Currently, the American Association for the Study of Liver Diseases, the European Association for the Study of Liver, and the Asian Pacific Association for the Study of the Liver recommend that HCC screening must be conducted at 6-month intervals [129]. But a consensus interval for recurrence after surgical resection has not been established. Moreover, this recommendation is based on the assumption that the HCC growth rate is similar in every patient. However, tumor growth in general is strongly affected by the microenvironment [130], and HCC growth rate is likely affected by host factors such as age, sex, preexisting diseases, etc., as well as by tumor factors, such as initial average HCC diameter, tumor multiplicity, etc. [131].
Therefore, we used our model to establish a proof of concept of how knowing different parameters that determine tumor development can predict such behavior. We have performed an exploratory sensitivity analysis by varying ±10% the input variables that feed our model: oxygen uptake of tumor cells and hepatocytes, hepatocytes and cancer cell cycles duration, repulsion and adhesion coefficients between cancer cells and hepatocytes. Figure 11 shows the tumor size relative change based on the variation of those parameters. The blue line represents the mean tumor size, the yellow shaded region represents the standard deviation of the mean value based on the stochasticity of the model. Red bars and blue bars represent tumor growth when the parameter original value is increased and decreased by 10% respectively. The parameters whose bars fall into the yellow zone do not modify the tumor growth, hence, according to our model only three parameters mostly control tumor growth. The largest influence is exerted by the hepatocyte oxygen uptake constant.
When we increase this constant by 10%, the hepatocytes need more oxygen to keep metabolically active. That causes a reduction in the oxygen concentration in the microenvironment which leads to a reduction in the cancer cell cycle entry. It shrinks the tumor final size by 66%. On the other hand, when the hepatocyte oxygen uptake constant is decreased by 10%, there is more oxygen available in the microenvironment, and the final tumor size will be increased by 200%. Reducing the same constant in cancer cells, also impacts in the final tumor size, however, not as significantly as in the previous instance. The second most influential parameter is the cancer cell cycle duration.
A 10% reduction leads to an increase of the tumor final size by 140%. Conversely, increasing this parameter by the same amount will shrink the tumor size a 47%. In clinical practice, those parameters can be associated to patient specific factors. It has been reported that cell division rates decrease with age [132,133], so that younger individuals have a higher fraction of dividing cells than older individuals. In fact, although the removal of up to 75% of the total liver volume is feasible in a young patient (≤ 40 years of age) with normal hepatic parenchyma, it is suggested that resection must be more conservative in elderly patients (e.g., ≥ 70 years of age) [68]. Other studies show that regenerating liver consumed an increased amount of oxygen, especially during DNA synthesis after hepatectomy [113,134], that is why dynamic assessment of preoperative exercise capacity may be a useful predictor of short-and longterm postoperative prognosis [135]. Even though at its current stage, our model can qualitatively and quantitatively capture characteristics of the processes of liver regeneration and hepatocellular carcinoma recurrence, the data from markers that assess the patient's liver overall health could be used to feed the algorithm with specific patient factors, and in this way the model would easily transform into a workbench for hypotheses testing.

Discussion
In this study, we developed a 3-D off-lattice agent-based model to simulate large multi-cellular systems with the ability to test tissue-scale effects resulting from cell-scale interactions. One of the main characteristics of this type of models is their high predictive power. Obviously, in this context, a pre-requisite is a careful calibration and validation of the ABM. This means, one should first tune the model with real data until the known tissue behavior is reproduced. Since the liver has a remarkable capacity to regenerate, liver regeneration after partial hepatectomy has long been an excellent testing ground for modeling tissue regeneration. Moreover, although the liver consists of various types of cells, hepatocytes account for about 80% of liver weight and about 70% of all liver cells, which simplifies some of physiological aspects the model has to address. One can focus only on hepatocytes in order to study the relation of organ size with number and size of cells, and this makes them the ideal candidates to be the agents of our model. On the other hand, taking into consideration that the liver has a double supply of oxygen through the portal vein and hepatic artery, intrahepatic blood flow is highly regulated both by the disposition of its anatomical unit, which is the lobule, as well as by the interaction of its components with the extracellular matrix.
Dirichlet nodes turn out to be the perfect computational construct to represent these blood vessels' architecture and function, and thus add further realism to the oxygen diffusion within the model organ.
It is to be stressed that the precise mechanisms of initiating, promoting and terminating regenerative responses remain to some extent unknown. Here we have proposed a computational regulatory mechanism based on a substrate that diffuses in the cell microenvironment, substrate here denoted as growth factor. We assume that when the liver is intact, the microenvironment is in homeostasis and, although they are metabolically active, hepatocytes remain dormant in the cell cycle. When a partial hepatectomy takes place, each cell secretes a constant amount of growth factor and, depending on its concentration, they will undergo either hypertrophy or proliferation, We consider important to stress the fact that despite the major simplifications that our ABM presents, its outcomes show that the overall functionality of the model is preserved. Since the liver does not recover its original shape, the reduced spherical model served properly to reproduce the liver regeneration dynamics. Moreover, since the grid we used to simulate the liver is orders of magnitude smaller than an actual liver, the region could be considered as representative of subregions of the liver, avoiding architectural issues. On the other hand, our results show that by tuning the effective oxygen diffusion and decay constants, we can reproduce the effect of sinusoids in the liver lobules. Although our model does not place limitations on the liver size and shape, and allow us to introduce some inhomogeneities in the microenvironment in order to simulate the sinusoids, we believe that it will increase the computational burden of the simulations and it will not make much of a difference to the results. The comparison between in silico and in vitro systems shows that even a model based on simple rules governing cell cycle, intercellular bonding and basic physical relationships between neighbouring cells can successfully reproduce the behaviour of a real biological system.
Once the model is calibrated, it can be used to study the emergent behavior of the tissue in different scenarios. First, we implemented our model to test another degree of PH. Since 50% PH has now been successfully and frequently used for surgical removal of liver associate tumors, we have performed an in silico 50% PH. Here, similar to the 70% PH, hepatocytes enlarge their volume and also proliferate but with a reduced replication rate. We adjusted the results of liver regeneration time and percentage of proliferation based on the degree of PH of our simulated PH, to get an approximation of the liver behavior under other resection degrees.
Second, since the PH is the first line of treatment for patients with hepatocellular carcinoma in stage 0-A of BCLC staging system without clinically significant portal hypertension, we used our model to predict the potential recurrence of the tumor in the remnant liver after a 70% PH. The outcome of our simulations is in accordance with clinical observations, which comes to reinforce our confidence in the applicability of this approach in other scenarios. In that context and based on the fact that HCC growth rate is likely affected by host factors, we have performed an exploratory sensitivity analysis by varying ±10% the input variables that feed our model. We found that there are three parameters that according to our model most likely make an impact on the tumor growth: the hepatocyte oxygen uptake constant, cancer cells cycle duration and cancer cells oxygen uptake constant. We take this as a proof of concept of how knowing different parameters that determine tumor development predict such behavior in this 3-D off-lattice agent-based model. We have also studied the tumor growth rate based on the initial position of the the residual tumor clone. We observed that tumor volume vary depending on the location of the residual tumor clone. It grows faster when is located at the periphery of the liver. This result suggest that, although the geometry of the liver does not influence when study liver regeneration, it might make an impact on the HCC recurrence.
It is important to mention that even though at its current stage, our model can qualitatively and quantitatively capture characteristics of the processes of liver regeneration and hepatocellular carcinoma recurrence, it is not calibrated to any particular type of patient specific parameters (age, sex, pre-existing conditions,...). This is obviously a handicap for the model's direct application in clinical practice. However, the algorithm can be fed with specific patient factors, and in this way the model would easily transform into a workbench for hypotheses testing. The data required for model calibration should include tissue biopsies, as well as data from markers that assess the patient's liver overall health. Such data set can be used to both quantify different cell types, and record their spatial arrangement. Half the images can be considered training data and be used to calibrate the ABM, while the other half can serve as testing data reserved to check if the model predicts tissue behavior with reasonable accuracy. In a way, this procedure resembles the training of a neural network.
Some limitations should be stressed in our current model. Firstly, the sample sizes are considerable smaller than those expected in real situations. As mentioned, in order to accommodate the simulations to temporal and spatial scales accessible to regular computer facilities, the lattice we used to model liver regeneration is 1 mm 3 in volume, far too small when compared with actual human liver size. Even if the region could be considered as representative of sub-regions of the liver, is obvious that important size effects might be left out. This challenge can be tackled resorting to new programming methodologies based on massively parallel computational approaches for multi-core and multi-tensor core devices such as the modern general purposes graphic processing units. Research along these lines is currently in progress.
A second limitation at the cellular-tissue scale is that some important cell types found in the liver, such as, biliary epithelial cells, stellate cells and Kupffer cells, are not included in the current model. Even though hepatocytes are responsible for most of the metabolic and synthetic functions of the liver, a future improvement could be the implementation of the interaction between different cell types and explore how it affects liver regeneration processes. A substantial amount of physiological research would be required as a prior step to the implementation of this model upgrade.
Finally, a few words concerning our model's implementation. We have been particularly careful to construct a software platform in a modular and extensible fashion. The aforementioned modules can be replaced with ones with more fine-grained versions as discussed, so that more specific details can be incorporated (as properties) and new processes as well (as methods) with different degrees of detail. Even though our model is not a 1 : 1 in silico copy of the liver and, therefore, it can not accurately describe in full detail the complex biology of liver regeneration and HCC recurrence, it could serve as a tool to test different hypotheses, as well as for testing and analyzing possible outcomes using multiple plausible parameter combinations. We are confident that once the goal of implementing patient specific factors is reached and the model undergoes a rigorous calibration and validation, it could be used as a platform for in silico conducting virtual clinical trials.