A multiscale agent-based model of ductal carcinoma in situ

Objective: we present a multiscale agent-based model of Ductal Carcinoma in Situ (DCIS) in order to gain a detailed understanding of the cell-scale population dynamics, phenotypic distributions, and the associated interplay of important molecular signaling pathways that are involved in DCIS ductal invasion into the duct cavity (a process we refer to as duct advance rate here). Methods: DCIS is modeled mathematically through a hybridized discrete cell-scale model and a continuum molecular scale model, which are explicitly linked through a bidirectional feedback mechanism. Results: we find that duct advance rates occur in two distinct phases, characterized by an early exponential population expansion, followed by a long-term steady linear phase of population expansion, a result that is consistent with other modeling work. We further found that the rates were influenced most strongly by endocrine and paracrine signaling intensity, as well as by the effects of cell density induced quiescence within the DCIS population. Conclusion: our model analysis identified a complex interplay between phenotypic diversity that may provide a tumor adaptation mechanism to overcome proliferation limiting conditions, allowing for dynamic shifts in phenotypic populations in response to variation in molecular signaling intensity. Further, sensitivity analysis determined DCIS axial advance rates and calcification rates were most sensitive to cell cycle time variation. Significance: this model may serve as a useful tool to study the cell-scale dynamics involved in DCIS initiation and intraductal invasion, and may provide insights into promising areas of future experimental research.


S1.2 Phenomenologically determining estimates of model parameters
In the calibration of our model, we encountered several mathematical quantities that we were unable to quantify from literature-reported values. These included diffusion constants for AREG, FGF, estrogen, as well as binding or release rates for these molecules, or the concentrations which upregulate proliferation. To circumvent this issue, we estimated diffusion coefficients for these molecules based on molecules with similar structures and molecular weights, and used normalized concentrations for these molecules within the simulations performed. Specifically, binding and/or release rates (λi, equations 3,4) of AREG, estrogen, and FGF were determined using standard parameter estimation techniques (basically, calibrated values which together minimized the discrepancy between simulated and reported DCIS duct advanced rates). Baseline signaling thresholds were estimated under the assumption that the highest literature-reported DCIS invasion rates [4][5][6] occur under conditions where cell proliferation is uninhibited by unavailability of signaling molecules involved in upregulating proliferation. Testing of the model under circumstances where proliferation was completely uninhibited by any signaling effects revealed that ductal invasion rates of DCIS are capable of exceeding those reported in the literature, and this thresholding plays an important part in regulating cell proliferation and associated duct advance rates. This allowed us to establish a baseline threshold for each molecule where literaturereported invasion rates were observed, and at the same time where all phenotypes were able to continue proliferation within the same simulation. We refer to this baseline signaling threshold as "medium" in Figures 7 and S2.

S1.3 mesh discretization, computational costs, and numerical details
Mesh discretization was on the order of the mature agent diameter, with element length around a median range of 10 µm. While it is not practical to enforce an exact element length across the entire mesh, the simple geometry of the duct domain allowed us to ensure mesh element length was consistent across the entire domain (as there was no need for any localized mesh refinement, etc.). Continuum solutions were obtained via the FEM method using direct linear solvers, 2nd order Gaussian Quadrature, L 2 norm, and 1st order Lagrange (scalar) basis at 60 second intervals, while ABM solutions were advanced in 30-minute intervals. Simulations were executed in serial (single node single core) runs on the Lonestar 5 computer at Texas Advanced Computing Center [7], with all results presented from simulations run for ≤ 48 wall clock hours (except parameter perturbation results in section 3.4, which were run up to 7 days wall clock time).

S2.1 Cell-cell physics
In our DCIS model, agents are bound to remain within the duct cavity (defined by the luminal inner layer of the mature duct), and may not transition to invasive ductal carcinoma via penetration of the duct basement membrane. Within the duct, agents may move freely as determined by the physics of cell-cell interactions (solved using Bullet Physics, an open-source rigid and soft body physics engine implemented in C++), and all agent coordinates are updated after every cell conformational and position change event (i.e., mitosis, movement, growth, lysis-induced swelling and cell rupture, etc.). Cells may be displaced by events in neighboring agents (e.g., movement, proliferation, growth), and are infinitely compliant to displacement. Subsequent to a cell proliferation or growth event, the movements of the surrounding cells are solved through a series of dynamically refined (variable) time steps to mechanical relaxation (achieved through both BulletPhysics adaptive internal time steps and through imposition of a fine time step resolution in the physics engine world), allowing the cell physics to return to a new equilibrium configuration before the next proliferation event or advancement to the next simulated time step. In this way, we are able to achieve a high level of stability in the cell-scale discrete model, and avoid any unforeseen behaviors that may result in instability. For example, a large number of simultaneous proliferation events (i.e., all proliferation events scheduled to occur within an ABM time step will be executed at the end of the time step) which is possible taking 30 minute time steps, where all proliferation events that would have occurred within that 30 minutes are executed upon advancement to the next time step) could potentially result in artificially large displacements (and thus large accelerations, which may incorrectly break cell-cell adhesion, etc.) in some cells due to the summation of many smaller simultaneous displacement events. By allowing a return to mechanical equilibrium after each event, we are able to avoid these potential complications and assure reliable stability in the discrete portion of the model without imposing any artificial constraints (for example, this can be accomplished artificially by restricting (clamping) movement to no more than a defined maximum distance per time step).
Cells adhere to neighboring cells through application of cell-cell adhesion forces, and are resistant to displacement due to these forces, as well as static, kinetic, and rolling friction (simulating cell adhesive forces). Cells may also deform due to pressure from surrounding cells, which we represent mathematically as a loss of energy through a coefficient of restitution CR, represented for each cell i as the square root of the ration of the cell's kinetic energy before and after deformation as where m is the cell mass and vc and uc are the cell velocity before and after cell-cell interaction, respectively (and thus cell deformations are not modeled explicitly as a change in plasma membrane conformation). However, because two cells may not occupy the same space, ultimately cells must be compliant to their neighbors, and thus displacement forces always overcome the forces which resist cell movement, although these forces do play a direct role in determining which displacements will occur.

S2.2 Maximum cell proliferation cycles
Normal mammary epithelial cells have been observed to divide 5560 times in culture, and even up to more than 250 cycles through immortalization with c-myc or its direct transcriptional target, hTERT (the human telomerase subunit responsible for catalysis) [8,9]. Although cancer cells are known to experience a level of immortalization, we choose to limit the maximum possible proliferation cycles (Pmax) of the DCIS progenitor phenotype to the lower end of this spectrum, i.e., Pmax ≤ 50 mitosis cycles before differentiation (to a terminally differentiated, i.e., nonproliferative state). DCIS progenitors count the number of cycles they have proliferated (akin to a telomere shortening mechanism), and upon hitting the proliferation threshold Pmax will always give rise to two terminally differentiated daughters. Cells are proliferation restricted by surrounding cell density, and will become quiescent if the local cell density (approximated by true mechanical density, i.e., cell mass/volume within a predetermined distance from each agent) crosses a certain threshold. This restriction is reversible, and if density of the surrounding cells is reduced due to cell death or movement, then a progenitor phenotype cell may resume proliferation. In the model, cells which have proliferated the maximum number of cycles will always give rise to two differentiated daughters in their mitosis cycle. We note that it is possible for Pmax to exceed the time simulated in the model runs presented here; in this case, the differentiated phenotype is not observed.

S2.3 Cell density-induced quiescence
We have implemented a somewhat simplified approximation for cell density as a density threshold over a certain number of "cell-lengths" (diameters) around the cell. There are many ways this can Taking into account the model simplification that cells are approximated as rigid-body spheres, we may use standard cell lattice sphere packing values (ignoring the fact that the sphere-packing becomes more complicated when spheres are not all equally sized) to obtain a maximum "actual" cell density, which is known to be 0.74 (FCC and/or HCP atomic packing factor) with equally sized spheres. Using this adjustment factor (i.e., that cell 100% cell density occurs at 74%), we may calculate the actual physical density (mass/volume) of cells within a certain number of cell diameters of the cell; that is, its "local density." Cell mass is determined by cell volume (for simplicity, we assume constant density across all cells), leaving only the question of what is the effective volume over which density should be calculated. That is, now many cell lengths away from the cell of interest does local cell density shown an effect? We estimated this through a calibration method against maximum observed DCIS advance rates measured in vivo; which we assume are representative of a case effectively uninhibited by the signaling pathways we examine. By testing the model over a range of these "effective density distances" under simulation conditions which allowed proliferation to be completely uninhibited by molecular signaling, we were able to determine a value that worked well to consistently replicate measured in-vivo observed DCIS advance rates.

S2.4 Hypoxia, necrosis, and cell lysis
Hypoxia is commonly seen in DCIS in ducts over a certain diameter; hypoxia-induced necrosis has been reported in up to 94% of ducts of diameter > 360 µm, but only 34% of ducts with diameter < 360µm, with an average viable rim thickness of 180 µm before hypoxia onset [4]. This is due, in part, to increased oxygen consumption rates in mammary cancer cells, reported to be 150260 Mol cell -1 s -1 for EMTGIRo and MCF-7 cell lines, respectively, relative to 2.545 Mol cell -1 s -1 for healthy cells (depending on cell type and cell phase, i.e., quiescence vs. mitosis phases) [3]; this translates into an increase of up to ~100 fold oxygen consumption in cancer cells. Cancer cells have been reported to enter hypoxia when local oxygen supply drops below 810 mmHg [10,11], and at about 1/3 the normoxia concentration observed in the healthy tissue [10], which is the threshold value we used for the onset of hypoxia in the work presented here (i.e., hypoxia threshold θH = 1/3 normoxia). Cells in the model will enter irreversible hypoxia induced necrosis after a programmed time to necrosis (τN = 12 hours of continuous hypoxia) based on values reported in [12]. Subsequent to necrotic death, cells enter a lysis phase, swell under lysis conditions until the plasma membrane ruptures. In the work presented here, time to rupture τL = 6 hours as per [13] (mammalian cells may swell from 15× up to at least 10× their original volume during lysis [14][15][16]; we use 2× in this work, as seen in previous DCIS modeling work [17]), resulting in their cytoplasmic contents being released into the luminal cavity. Leaked cytoplasmic contents may then form microcalcifications (as seen in mammographic imaging); here we take the calcified volume (Vc ) to be 30% of the volume of the cell at the time of lysis initiation (as cytoplasm is reported to be roughly 70% water [18]). Spilled cytoplasmic contents become calcified after a programmed time to calcification τc =14 days (calcification times in 4T1 and 4T1.B cell lines in vitro are reported as 11 and 14 days, respectively, with 4T1.B cells more commonly forming microcalcifications in vivo [19]; thus we choose 14 days for this study).

S3.1 Number of TIC cells in TIC niche
We assumed that the number of cells undergoing a transition to a TIC phenotype will be minimal, and may occur in one or a small cluster of adjacent cells (as lineage tracing studies of the mature mammary gland homeostasis maintenance process have shown that over time, many contiguous cells may share a common mother, clonal lineage, and genetic makeup [20]). Under this assumption, we tested the effects of initiation of a TIC phenotype in between 19 luminal cells (specifically, 1, 3, 5, 7, and 9 TICs) at t = 0 on DCIS population expansion rates and early duct advance rates, as shown in Figure 4 (main text).
By examining the extent of DCIS advance over time, we discovered that the number of TICs initiated in the niche at time t = 0 had a negligible influence on the total DCIS extent and rate of advance for all three duct sizes tested (i.e., 100, 150, and 200 µm in diameter). Average growth rates over the time were estimated by linear fitting of data shown in Figure 4, and were found to increase from 11.84 to 11.87 mm/year for the 100 µm duct, 12.82 to 12.89 mm/year for the 150 µm duct, and 10.78 to 11.21 mm/year for the 200 µm duct (a 0.25%, 0.55%, and 3.99% difference, respectively). This result occurs because the disease rapidly transitions from an exponential growth phase to a linear cell-density restricted growth phase, and demonstrates only minimal influence of cell number variations in the TIC niche on duct axial advance rates.

S3.2 DCIS axial advance rates are consistent over all duct sizes tested
To test the influence of duct size on the ductal advance rate of DCIS, we performed simulations in all three duct diameters tested under baseline values (3 simulations for each case), but without the effects of signaling thresholds, as shown in Figure S1. Because of the complex interplay between the distribution of cell phenotypes and the associated upstream (e.g., estrogen) and downstream (e.g., AREG, FGF) signaling (also see Figure 1, main text), DCIS advance rates under signaling threshold effects may experience significant variation. Thus, we remove this effect to eliminate interaction of these variables to better examine only the effects of variation of the duct diameter.
Interestingly, we observed that the variation of DCIS axial advance rates was small between the different duct diameters, with average ducal advance rates (as estimated by linear fitting) averaging from 11.23 mm/year for the 200 µm duct to 12.95 mm/year for the 150 µm duct (and a more centralized rate of 11.91 mm/year for the 100 µm duct). This relatively small variation of axial advance rate relative to the variation in duct diameter is largely due to the effects of cell density in this case, and is examined more thoroughly in Discussion, main text.

S3.3 Further examination of signaling-limited DCIS invasion
In order to further elucidate the effects of signaling thresholds on DCIS ductal invasion rates, simulations were repeated in triplicate at each of the high, medium, and low thresholds (i.e., Figure   8, main text); this is shown in Figure S2. It was observed that when thresholds were low ( Figure   S2, A-C), cell proliferation was able to proceed mostly uninhibited due to sufficient endocrine and paracrine molecular concentrations throughout the simulated duct regions. When signaling thresholds were high, proliferation inhibition within the DCIS leading edge due to insufficient signaling stimulation was found to be dependent on the phenotypic makeup of the leading edge.
In this case, the ability of the leading edge to continue to advance was dependent on the phenotypic makeup of the leading edge, where presence of the uninhibited phenotype allowed for continued proliferation within that sub-population; this was manifest as slower total axial invasion and larger standard deviation among simulation runs (Figure S2, G-I). Moreover, it was observed that restriction in the estrogen signaling caused greater total reduction in axial advance rates, as estrogen acts upstream of the AREG-FGF pathway (this is discussed further in main text). At the medium thresholding level, significant proliferation inhibition was not observed in the simulations performed ( Figure S2, D-F).  as mean total axial ductal invasion distance, with standard deviation across all runs indicated by error bars. High thresholds were observed to cause significant inhibition of axial ductal invasion rates, with estrogen-limited profusion having greater effects than FGF-limited diffusion. This is likely due to estrogen acting upstream of FGF, and is explored further in main text.