Multicompartment modeling of protein shedding kinetics during vascularized tumor growth

Identification of protein biomarkers for cancer diagnosis and prognosis remains a critical unmet clinical need. A major reason is that the dynamic relationship between proliferating and necrotic cell populations during vascularized tumor growth, and the associated extra- and intra-cellular protein outflux from these populations into blood circulation remains poorly understood. Complementary to experimental efforts, mathematical approaches have been employed to effectively simulate the kinetics of detectable surface proteins (e.g., CA-125) shed into the bloodstream. However, existing models can be difficult to tune and may be unable to capture the dynamics of non-extracellular proteins, such as those shed from necrotic and apoptosing cells. The models may also fail to account for intra-tumoral spatial and microenvironmental heterogeneity. We present a new multi-compartment model to simulate heterogeneously vascularized growing tumors and the corresponding protein outflux. Model parameters can be tuned from histology data, including relative vascular volume, mean vessel diameter, and distance from vasculature to necrotic tissue. The model enables evaluating the difference in shedding rates between extra- and non-extracellular proteins from viable and necrosing cells as a function of heterogeneous vascularization. Simulation results indicate that under certain conditions it is possible for non-extracellular proteins to have superior outflux relative to extracellular proteins. This work contributes towards the goal of cancer biomarker identification by enabling simulation of protein shedding kinetics based on tumor tissue-specific characteristics. Ultimately, we anticipate that models like the one introduced herein will enable examining origins and circulating dynamics of candidate biomarkers, thus facilitating marker selection for validation studies.

Although the search for blood-borne cancer protein biomarkers to improve diagnosis and prognosis has accelerated in recent years, only a handful of candidates have reached clinical application [1][2][3][4][5][6][7] . A major reason is that the relationship between protein abundance in tumor tissue and in blood circulation remains poorly characterized. It has proven difficult to link the cell-scale events (e.g., protein shedding into circulation) to the dynamically evolving tissue-scale conditions (e.g., tissue access to vasculature that is changing as a result of proliferating and necrosing tissue). As the cell activity (proliferation and death/necrosis) and the associated protein shedding occur on similar time scales, the detection of proteins in circulation presents primarily a spatial problem, as proteins must persist and diffuse through space from cells to blood under conditions that are continuously changing this space and the access to the blood.
To complement experimental efforts, mathematical modeling and computational simulation techniques have been applied to predict circulating biomarker levels from tissue-scale data [8][9][10][11][12][13] . Previous continuum-scale modeling of tumor growth has considered the effects of molecular and cellular heterogeneity on tumor spatial progression (e.g., 14,15 ), but not in the context of blood-borne biomarkers. In particular, Gambhir and coworkers developed a 1-compartment model to simulate shedding of the secreted ovarian cancer marker CA-125 into blood plasma 8 . Applying this compartmental approach and its assumption of tumor and healthy cell populations contributing uniformly to shedding, tumor tissue and its vasculature were modeled as singular homogeneous Scientific RepoRtS | (2020) 10:16709 | https://doi.org/10.1038/s41598-020-73866-8 www.nature.com/scientificreports/ entities without cellular level detail. The strength of this model lies in a small set of parameters that can be easily interpreted and tuned. However, the model over-simplifies the tumor growth process, neglecting intra-tumor heterogeneity due to inhomogeneous vascularization that directly affects marker shedding into the tumor microenvironment and the blood circulation. Additionally, the model lacks the ability to simulate non-extracellular (non-EC) proteins (surface or intracellular in domain). As less than 20% of cellular proteins are secreted 16 , it may be challenging to explore the majority of candidate biomarkers using this approach. Building upon the work of Gambhir and coworkers 8 , we present a new model that links the levels of protein shedding into circulation to the heterogeneous tumor spatial vascularization and the corresponding proliferating and necrosing cell populations. This is achieved by discretizing tumor tissue into dynamically evolving compartments based on cellular distance from vasculature 17 , effectively bridging from the cellular to the tissue scale, and thus enabling prediction of system-level effects from cellular-scale events. In this model, proteins shed by cells diffuse through interstitial space to reach the vasculature in order to enter the blood circulation. It is well known that incipient tumors begin to slow growth and require sustained angiogenesis to meet the demands of their growing cell populations 18 . Vascularization induces gradients of oxgygen and nutrients in surrounding tissue, imposing cellular stresses such as hypoxia-induced necrosis, which influence shedding of both extracellular and non-extracellular proteins 9 . Intratumoral heterogeneity is induced by cells with higher net growth rates in well oxygenated regions proximal to vasculature, and with lower net growth rates in hypoxic and necrotic regions distal to the blood supply 19 . By simulating the dynamics of tumor development and the spatial variation in cell viability across heterogeneously vascularized tumor tissue, the shedding kinetics of a wide variety of proteins can be described as a function of tumor tissue characteristics. In particular, by enabling an explicit quantification of necrosis, it becomes tractable to model the shedding of the tens of thousands of non-extracellular proteins.

Results
A compartmental model to represent protein shedding by heterogeneously oxygenated tumor tissue. This study focuses on two fundamental coupled features of tumor tissue that impact protein shedding kinetics: cells that proliferate or die (necrose), and heterogeneous vascularization. We design a model to simulate proliferating and necrosing populations. Modeling these two populations distinctly allows for differential shedding rates of extracellular proteins and non-extracellular proteins. Specifically, viable cells dominantly release extracellular proteins, whereas necrosing cells dominantly release intracellular proteins through cell membrane decay and decreased production and secretion of extracellular proteins. Once exported from the cell, these proteins must diffuse through interstitial space to reach the vasculature. Heterogeneous vascularization leads to cell subpopulations experiencing varying levels of oxygen and nutrient availability (Fig. 1A). Differential oxygen and nutrient availability in turn impact cellular proliferation and necrosis, and the associated release of proteins into the interstitium. This allows modeling the difference in shedding rates between extraand non-extracellular proteins from viable and necrosing cells as a function of heterogeneous vascularization.
To achieve this dynamic representation of protein shedding kinetics, the tumor model is constructed as a series of compartments. At any one time, we assume that each compartment acts as a source of proteins and as a sink for oxygen and nutrients. The whole vasculature within the tumor is represented as a source of oxygen and nutrients and as a sink for proteins. We assume that tumor cells can be grouped into distinct compartments in terms of their distance to vasculature, with each compartment having similar proliferating and necrosing characteristics. Thus, the model virtually rearranges tumor tissue into a set of compartments as a function of access to the vasculature (Fig. 1B). The vasculature is represented as a cylindrical (tube) compartment, surrounded by concentric cylindrical (shell) compartments of tumor tissue. Each compartment represents the subpopulation of tumor cells being at any one time at the same distance from vasculature. Inner compartments represent tissue proximal to vasculature, while outer compartments simulate tissue distal to vasculature. Tumor and vascular growth are captured by longitudinally extending these compartments along the cylindrical axis as a function of time.
In this manner, vascular-induced intratumoral heterogeneity is represented in the model as a function of radial distance to proximal vasculature, which determines cellular proliferating and death rates in differing oxygenation conditions. It is well known that tumor cells are more proliferative in highly oxygenated conditions and more necrotic in hypoxic conditions, e.g., as shown by cell line data collected from non-small-cell lung cancer cells in vitro, which we use here to set compartmental cell birth ( k B ) and death rates ( k D ) 20 (Fig. 2). In turn, the cellular proliferation and necrosis in each compartment affects the respective shedding of extra-cellular and non-extracellular proteins. The proteins are shed into the centrally located vascular compartment according to a weight (w) representing diffusivity based on the compartment's distance to the vasculature (Fig. 1C). With this compartmentalized representation of tumor vascular-induced heterogeneity and its consequences on protein shedding, the model can be formulated with constraints based on observable tumor characteristics. In particular, tissue spatial constraints, e.g., as can be observed from histology 21,22 , are used here as model parameters.
The key parameters that link the spatial relationships in the simulated tumor are chosen as the following ( Fig. 1): Mean vessel diameter: The mean vessel diameter (MVD; denoted as d) is the average diameter of all vascular structures (e.g., blood vessels, capillaries, microvessels, etc.) observed across tumors. The average diameter of larger vessels that sprout microvessels and capillaries has been measured to be approximately 60 µ m across a set of histology slides of retinoblastoma 23 , which is used here as a representative value. The MVD varies between tumors of different tissues.
Relative vascular volume: The relative vascular volume (RVV) is the ratio of vasculature volume to tumoral volume. This metric has been observed to remain steady during tumor development, e.g., a value of approximately 17% 24 has been measured in mammary adenocarcinoma 25,26 . This value can vary based on the tumor tissue type. Motivation behind modeling heterogeneous subpopulations in the tumor microenvironment. The model approximates tumor growth as vascularization occurs early in development (abstracting various methods of tumor vascularization, including angiogenesis, vasculogenic mimicry, and microvessel formation). The microenvironment cross-section represents any given 2D neighborhood perpendicular to the vasculature. The vessel diameter d, necrotic cuff ǫ , and cylindrical model radius R are parameter values derived from cellular spatial constraints obtained from histology images 23 . The radius R of the cylindrical model is equal to the combined length of half the mean vessel diameter ( δ = d/2 ) and the necrotic cuff ( ǫ ). Thus, the cells in the tumor are approximated by the sum of all such 2D cross-sections taken over the length of the total vasculature as governed by h t . Tumor cell subpopulations are defined by their radial distance to proximal vasculature due to heterogeneous access to oxygen. This radial parameterization is discretized with compartments that correspond to each cell subpopulation. Tumor compartments vary in volumes and carrying capacities and are calculated with the recurrence relation of subvolumes from cylindrical shell and spatial constraints described herein. The kinetics of vascularization updates the constraints on compartmental volumes and carrying capacities at each time step. (B) Compartment diagram and parameters for tumor growth. Parameter σ is the uniformly partitioned thickness of each compartment, set at approximate single-cell diameter of 10 µ m. Compartments T i all generate cells with their specified birth, death, and net growth rates ( k G,i ) based on their access to oxygen (i.e., distance to vasculature). These compartments thus experience proliferation and necrosis at differing rates, which subsequently leads to differing shedding behaviors. Cell motility with a preference toward oxygenated regions allows for added dynamism of the model (denoted as probabilistic terms between compartments). (C) Compartment diagram for protein shedding from tumor cells to plasma. Extracellular protein shedding is dependent on the net proliferative population over the compartments and (per-cell active shedding rate), while non-extracellular (intracellular and surface) protein shedding is dependent on the instantaneous cell death over the compartments and (per-cell instantaneous contribution due to turnover). Shedding to plasma, denoted as the Pl compartment, is further weighted by distance to vasculature with weights w i to account for heterogenous diffusion based on distance.
The model implements a longitudinally growing cylindrical tumor tissue of height h around a growing vessel. Here, ǫ is the radius of the cross-sectional model and the MVD constrains the vascular diameter. As tumor vascularization progresses in response to a net release of pro-angiogenic stimuli by tumor tissue, the tissue carrying capacity and growth correspondingly increase. A near-constant RVV allows for the assumption that the MVD and ǫ of the microenvironment cross-section remain constant as the carrying capacity of the tumor subpopulations increases by increased vascularization.
In this framework, a growing vascularized tumor is represented by n discrete tumor compartments T i of single-cell thickness that shed proteins into the vascular compartment Pl (denoting plasma) based on their tumor populations (proliferating or necrotic) and their distance to the vascular compartment.
Shedding Parameters: A key component of this model is its ability to capture the shedding of both extracellular and non-extracellular (cell membrane or intracellular) proteins. Prior models exclusively focused on extracellular proteins. However, experimental studies have clearly demonstrated that intracellular components are able to shed to the circulation. This model enables capturing the shedding kinetics of extra-and intra-cellular proteins from proliferating and necrotic cells in heterogeneously vascularized tumor tissue. Key parameters describing the shedding kinetics include , the protein-specific contribution per cell during necrosis, and , the proteinspecific shedding rate per cell per day. These two parameters operate in opposition to best capture the differences in shedding dynamics of extracellular vs non-extracellular parameters. Specifically, extracellular proteins will have low values of and higher values of that relate to secretion rate. Non-extracellular proteins will have higher values of calibrated in an abundance-dependent manner, and low values of . Variable u i (t) is defined as a distance-based relationship between the number of tumor cells and the corresponding compartment's population and shedding values to capture the diffusion of molecules from varied parts of the tumor into the blood. Variable q(t) is the final summarized extracellular and non-extracellular protein mass in plasma at time t.
Model calibration. The first step to employ the model is to choose a set of parameter values that yield tumor growth that is biologically realistic. For this purpose, we chose to match the tumor cell proliferation shown by Hori et al. 's model, which has been calibrated from experimental data with ovarian cancer cells 8 . To run the model over a range of possible parameter values, scans over a parameter set of both the maximum volumetric concentration of oxygen at the vasculature point ( C 0 ) and the vascularization rate ( k V ) were simulated. Parameter C 0 helps to determine the range of values for the n = 10 compartmental net growth rates k G,i , helping to shift the range of possible proliferation rates depending on the proliferative capacity of the particular cell type. Because of the large number of parameter combinations, validity cutoffs were arbitrarily defined (e.g., tumor cell proliferation in all n = 10 compartments, total cell population of greater than 10 7 cells at t := t end = 12 years) to identify suitable tumor growth simulations. These scanned parameter pairs are presented in Fig. 3. Tumor cell proliferation matching that of Hori et al. 's model was achieved given the model's linearity and spatial constraints, showing it to be a reasonable spatial extension of their dimensionless model. Using the subset of suitable results, the simulation with maximum growth was selected to evaluate protein shedding profiles.
Sensitivity of tumor growth and protein shedding to variation in parameter values. Next, the sensitivity of the calibrated model to variation in its parameter values was evaluated. Running these analyses over the full set of parameters introduced by the model, the influence of each parameter on tumor growth ( Fig. 4) and protein shedding (Fig. 5) was evaluated for both non-extracellular (non-EC) and extracellular (EC) proteins. These studies examined how stable the model was to variations in or mis-estimations of parameters. For tumor growth, great sensitivity is seen to maximum oxygenation concentration C 0 and vascularization rate Tumor compartment birth and death rates. Extrapolated rates are derived from cells exposed to varying levels of hypoxic stress 20 . The net growth rates ( k G,i ) are equal to the birth rates ( k B,i ) minus the death rates ( k D,i ). A nonlinear relationship can be observed as a function of compartment number (related to distance from proximal vasculature).

Scientific RepoRtS
| (2020) 10:16709 | https://doi.org/10.1038/s41598-020-73866-8 www.nature.com/scientificreports/ k V . Parameter C 0 modulates the slope of total tumor growth since the parameter set of local net growth rates is shifted by the interpolation function as described in Methods. As a result, constituent compartments are assigned local net growth rates that are shifted accordingly. Parameter k V appears to enforce a ceiling function for the total tumor population after it outpaces vascularization and necrosis occurs at higher rates. For protein shedding, t 1/2 (protein half life in blood), C 0 , and or all carry early influence on model performance by large orders of magnitude given realistic ranges specific to each parameter. One may consider both k V and C 0 hyperparameters tuned to a particular tumor, with considerations such as cell and tissue type, while or and t 1/2 are parameters that can be estimated a priori for a shedding protein of interest. For these analyses, the performance of simulated non-extracellular and extracellular proteins was compared. Non-extracellular trajectories tended to be lower and start later than extracellular protein trajectories. This difference stems from the source of non-extracellular proteins being dominated by necrotic cells. Accordingly, the This surface utilizes both of the parameters k V (vascularization rate) and C 0 (maximum volumetric concentration of oxygen in the tumor, located at vasculature) to form a coordinate mesh along the x-and y-axes. The z-axis represents the final proliferative population at the userdefined maximum iteration for simulation, t end . Each tile on the surface represents a simulation with the corresponding parameters on the x-and y-axes. If the value of k V is large and C 0 is too small, cells primarily stay in compartments near vasculature due to large carrying capacities and little outgrowth. If the value of k V is too small and C 0 is large, carrying capacity is not increased quickly enough and tumor cells die off too quickly to reach distant compartments. Accordingly, valid tumors are defined as those supporting heterogeneous subpopulations and a large overall population. Arbitrary cutoffs were defined to help select for such tumors: those that experience growth in all n = 10 compartments (indicated by the color bar) and grow to a population of at least 10 7 cells (represented by translucent plane) at the set maximum iteration of t end = 12 years. Tumors used for downstream simulation were both above the 10 7 population plane and with n = 10 populated compartments, but were also chosen for having the smallest population sizes given the aforementioned selection criteria to select for more necrotic tumors. A sample of five trajectories were taken to generate the lines seen in panels ( www.nature.com/scientificreports/ relatively small size of the necrotic population, as compared to that of the proliferative population (Fig. 3), early in tumor development leads to a lag in the shedding of non-extracellular proteins. Despite this lag, our analysis showing overlapping trajectories confirms that non-extracellular proteins can be accessible in the circulation and may reach abundances that exceed those of some extracellular proteins.
protein shedding as a function of tumor cell proliferation and necrosis. Next, we evaluated the capability of the model to simulate tumor proliferative and necrotic population evolution in time and the corresponding protein outflux. A goal of the proposed model is to provide a mechanistic understanding of the shedding kinetics of tumor cell subpopulations. Figure 6 shows the per-compartment populations of both proliferative and necrotic cells. In panel (B), an inflection at approximately six years of simulation time is noted. This is due to necrotic cell populations outpacing the rate of vascularization. Spatially, and in the context of the presented model, cells begin to occupy more necrotic compartments (distal to vasculature) due to the linear vascularization function. Figure 7 shows a similar effect for the cell subpopulations' marker outflux and contributions to overall shedding. Changes in rate for both panel (A) and (B) can be seen at approximately 10 5 and 10 2 cells respectively, indicating that the simulated tumor reaches a steady-state in marker shedding at later stages.

Model validation.
As an initial validation, the model-simulated protein shedding was compared to that of Hori et al. 's experimentally tuned model (Fig. 8). The results show that the proposed model can be used as a spatial extension of one-compartment models, indicating similar trajectories between the proposed model and the Hori shedding model for an extracellular marker such as CA-125. However, there is a notable difference in late-simulation behavior due to the linear vascularization function growing too slowly for the total tumor population. This slowdown causes increased necrosis and increased non-extracellular shedding, which is a model-derived hypothesis to be tested in vivo. With the selected extracellular protein features from CA-125, hypothetical non-extracellular protein shedding trajectories were generated. The comparison against this extracellular marker was performed via parameter scanning over possible values of protein-specific parameters , u H (healthy cell influx rate of protein mass), and t 1/2 . The colored regions show the dynamic range of trajectories for each value of u H , which effectively sets the basal mass that the growing tumor builds upon. The dynamic range is determined by the other two protein-specific parameters controlling plasma influx ( ) and elimination (as determined by t 1/2 ). Overall, this analysis indicates that non-extracellular proteins may outperform an extracellular protein for blood-based detection given boosted parameter values, highlighting the need for increased study of non-extracellular proteins for cancer biomarker discovery.

Discussion
This study implements a novel spatial modeling approach that enables in silico simulation of extracellular and non-extracellular protein shedding by proliferating and necrosing cell populations in vascularized tumors. The model may be useful for describing the origins of non-extracellular biomarkers in the blood, and creating a scalable framework for asking questions about the impact of environmental heterogeneity on biomarker shedding. These questions are critical, as accessing non-extracellular biomarkers opens up the possibility for thousands of additional proteins to be potential biomarkers. In addition, prior to the development of this model, to the best www.nature.com/scientificreports/ of our knowledge there has been no straightforward way to link variation in tumor vascularization to marker shedding. A key potential utility of the model will be in comparing the behavior of different candidate markers to each other for the purposes of selecting which markers to carry forward into validation studies. While the identification of such putative markers is beyond the scope of this study, we demonstrate in Fig. 8, that there could be non-extracellular proteins whose circulating abundance can exceed that of extracellular proteins. Individual proteins represent trajectories within these windows that can be directly compared.
Hori & Gambhir 8 proposed a general model to evaluate biomarkers shed by tumor and host tissue, without distinction as to their intra-or extra-cellular origin, in order to determine the limits of tumor detection from biomarker levels in blood circulation. The model does not consider tumor-specific characteristics, such as vascularization parameters, which can greatly influence the amount and type of biomarker outflux into circulation. Similar to Hori & Gambhir's model, this study employs a mono-exponential function to simulate tumor growth. Unlike their study, here we are interested in evaluating both non-extracellular and extracellular protein shedding kinetics during this growth, as a function of the tumor tissue characteristics. The model can be tuned by evaluation of histology data, including relative vascular volume, mean vessel diameter, and distance from vasculature to necrotic tissue. This information is used to simulate the dynamic interaction between tumor tissue and vascularization, which determines biomarker outflux into circulation.
The model discretizes cell populations sharing similar oxygenation as a function of distance to vasculature into interconnected compartments. To bypass the spatiotemporal complexities of vascularization, these compartments are in the model assigned increased carrying capacities as the tumor grows. The spatial representation is achieved through the use of concentric cylindrical compartments and, correspondingly, use of the cylindrical parameterization of radial distance to a central vascular compartment. Due to this multi-compartment approach, approximate solutions to difference equations are utilized because (i) smaller granularity than at the single day level is not needed, (ii) the time-scales of interest are large, and (iii) cellular proliferation in the tumor microenvironment is highly dynamic. For a given tissue type and its corresponding growth rates, the model can compute  Table 1. This analysis emphasizes the importance of all parameters in shedding of both EC and non-EC makers, as seen by their visibly large log-fold changes in protein mass. Namely, parameter u H controls the immediate uptick in the trajectory, while parameters t 1/2 and or control the slope of the trajectory after proliferative growth slowdown and the emergence of the necrotic population. Parameter t 1/2 also controls the initial surge in marker mass before reaching steady state due to controlling the rate of elimination. It should be noted that varying either k V and C 0 alone in the specified parameter ranges does not visibly benefit non-EC shedding since cell death is the source of markers. www.nature.com/scientificreports/  Consequently, this approach enables the comparison of the shedding kinetics of multiple proteins at once given a tumor simulation calibrated to particular tissue characteristics. The implementation is computationally low-cost, providing for spatial modeling without the need to solve coupled partial differential equations 21 or to represent individual cellular agents 28 . Due to the model's assumption of coupling between tumor (proliferative and dead) cell populations and plasma protein concentration, as well as an assumed lack of feedback from the plasma and tumor compartments, the tumor growth and protein shedding systems of equations can be computed in series. This computationally efficient methodology leverages local behavior of intratumoral dynamics to simulate the protein shedding dynamics of whole tumors. Accordingly, the proposed model offers a novel framework for rapid simulations of tumor shedding with large-scale comparisons between proteins of interest, with the potential to provide insight into patient-specific tumor conditions. For simplicity, the model employs a standard mono-exponential growth function to describe compartmental growth, yielding relatively small parameter sets both when generating estimators and when calibrating the model. Other growth functions could be explored. Moreover, dying cells (e.g., via necrosis or apoptosis) at each time step do not take up physical space towards the capacity of the compartment and, instead, are instantaneously cleared out of the system. Future tuned parameter sets of net growth rates will need to account for this assumption of space availability. Further, the model parameterization scheme presents certain limitations. The model assumes that cell are sourced from oxygenated compartments at time t = 0 , but in reality there is some a priori structure for populations before new vasculature enters the system through angiogenesis or vascular recruitment. Additionally, cell death is not solely caused by necrosis; while increased death rates in less oxygenated compartments imply hypoxia-induced necrosis, the proportion of cell death caused by necrosis as opposed to scheduled cell death is unknown. Furthermore, the vascularization function is currently assumed linear with respect to time. While a density dependent function would reflect angiogenic factors having a role in angiogenesis, a simpler model was chosen. Figure 8   www.nature.com/scientificreports/ function controls the rate at which cells can proliferate due to increased carrying capacity. Tuning is required to ensure that the function is not too fast (i.e., large-scale necrosis is delayed) nor too slow (i.e., the tumor fails to achieve a large enough population over time). Values for the model parameters were primarily chosen from the literature (Table 1). In particular, net tumoral growth, birth, and death rates were tuned to previously-obtained experimental data by us and by others 18 , serving as order of magnitude approximations. The MVD (and resulting volume of vasculature) was chosen as an approximate baseline. Additionally, k V was itself tuned to achieve a similar growth trajectory as presented in Hori et al. These tumor-specific values should be tuned to a particular tissue type of interest. In order to further tune and validate the proposed model, future work will aim to grow in vivo xenografts composed of cell types that shed selected proteins of interest. Known protein tumor kinetics will need to be used to validate model predictions and to verify its utility in order to bring it closer to clinical application. With repeated blood sampling and serial measurements (similar to studies such as Li et al. 29 ), the shared tissue-specific model parameters (e.g., growth and vascularization rates) can be tuned to the desired tissue type. After such tuning, protein-specific parameters ( , , u H , and t 1/2 ) can either be estimated using approaches that average over sampling time or calculated via mass-balance equations once xenografts have nearly reached carrying capacity and steady state is assumed. Lastly, future work will seek to build statistical learning models to predict protein-specific parameters 30,31 from proteomic features, publicly-available datasets 32 , and steady state values in order to fully enable large-scale viability studies and prioritization of candidate biomarkers. where, P(t) is proliferative population at time t ∈ R ≥0 and 0 ≤ k G ≤ 1 is net tumoral growth rate (Table 1). For model flexibility and cellular dynamics at each time step, difference equations are used at t: = 1 day intervals  (Fig. 3). & Values/ranges simulated for tumor growth sensitivity analyses (Fig. 4). + Values/ranges simulated for marker shedding sensitivity analyses (Fig. 5). # Values/ranges simulated for parameter scan comparisons between EC and non-EC shedding (Fig. 8). $ Values/ ranges simulated for CA-125 shedding simulations (Fig. 8). www.nature.com/scientificreports/ to form the Malthusian expression, P t+1 = P t + P . Through this use of a sufficiently small t , Eq. (1)'s differential equation is converted to Eq. (2) via a forward Euler finite difference approximation:

Methods
where P t is the net proliferative population at discrete time point t ∈ N and t = 1 day is the time step taken for this approximation. The initial condition is set to P t := 1 to signify the number of tumor cells that initiate tumorigenesis. In order to approximate cell death at time t, the net tumoral growth rate k G is decomposed into its constituent birth rate k B and death rate k D : Using the above definition of net growth rate, the terms in Eq. (2) are expanded to define B t and D t , the amount of cell birth and death at time t, respectively: where D t is the instantaneous cell death at discrete time t, encapsulating both regularly-timed cell death (i.e., apoptosis) and death by injury (e.g., hypoxia-induced necrosis). However for simplicity, D t is used an approximation of intratumoral necrosis, since this phenomenon governs much of the cell death in hypoxic regions 33 . This formulation for mono-exponential tumor growth and corresponding cell death will be used to describe the growth kinetics of each compartment's subpopulation, as described below.
tumor compartmentalization by radial distance to vasculature. We introduce a cross-sectional representation of the tumor microenvironment into compartments of shared tumor growth based on radial distance to proximal vasculature. With this radial dimension, a cylindrical model extends the cross-section to represent the growing tumor, as discussed later. The oxygen gradient maintained radially from the vasculature is parameterized and discretized to compute the numerical solution. The discretized regions define tumor cell populations with the same proliferative and dying potentials, defining n tumor cell compartments T i , each with its own growth rate k G,i as a function of distance to vasculature (Fig. 1). These compartments form a coupled system of linear ordinary differential equations. It should be noted that parameter σ , or the partitioning resolution, is the compartment thickness of the radially-partitioned oxygen gradient. To model subpopulations at the highest resolution, σ is set to approximate a single-cell diameter at 10 µ m. Along with a necrotic cuff of 100 µ m, this discretization yields n = 10 distinct compartments or tumor regions. The following formulation is used to define the radial relationship between compartments and their growth rates:

Definition 1
The level set corresponding to compartment T i of a real-valued function f of m variables is defined as a set of the form: where i ∈ Z.
In order to parameterize the system, points in R 3 are mapped to cylindrical coordinates, giving rise to the spatial dimension of radial distance to proximal vasculature, r. Accordingly, the net tumoral growth rate k G is defined as a function of r: where the coordinate (x, y, z) is a point in R 3 and r x 2 + y 2 is a point along the now parameterized radius. By using the above parameterization, Definition (1), and the assumption that cells will share a similar growth pattern locally, the system is descretized by level sets that define distance-based oxygenation regions of the microenvironment cross-section's constituent tumor regions and central vessel: where ⌈·⌉ is the ceiling function and f is a discretizing function that maps cells of a radial distance to a discrete region, or concentric shell of the tumor. As discussed in more detail in the following sections, the ith tumor cell compartment ( T i ) is assigned a level set that defines its corresponding volume v i and net growth rate k G,i . Together, these parameters encode for the distinct tumor regions to form a multi-compartment system of tumor growth. With this discretization, the net growth rate is defined as a discrete function rather than its continuous counterpart in Eq. (5). Thus, parameter k G,i is referred to as the local net growth rate: This mapping is used to define compartmental rates and shedding weights based on the corresponding radial midpoint.
compartmental volumes. The radial parameterization and discretization of space in the microenvironment into compartments also requires a spatial definition to encode heterogeneous population growth and shedding. With a radial dimension, a cylindrical representation defined around vasculature is used. The previous subsection introduced the level set formulation to assign compartmental rates to each subpopulation based on distance to proximal vasculature; this subsection introduces geometric definitions of cylindrical shells to define compartmental volumes. These two ideas together delineate constraints for both the tumor cell subpopulations and the vasculature (denoted as Pl in Fig. 1).

Definition 2
The volume V of a cylinder with radius R is:

Definition 3
The core C of a cylinder is defined as the sub-cylindrical region starting at the origin with radius δ.

Definition 4
The ith concentric shell (corresponding to compartment T i ) of a cylinder has a volume V i of: Definition 5 A cylinder with radius R and core C is said to contain n equidistantly concentric shells if its radius is uniformly partitioned by σ = R n . Namely, σ is the thickness of each cylindrical shell (i.e., associated with compartments T i ∀i = 1 . . . n ), where the ith concentric shell has radius: Furthermore, the shell subvolumes of v i (∀i = 1 . . . n) s.t. n i=1 v i = V n = V are computed in the following manner: Proposition 1 For a cylinder with radius R and equidistantly concentric shells, the following recurrence relation and solution describe corresponding subvolumes v i : Given that the MVD (d) has been measured at approximately 60 µ m in some tissue types, the radius of the average vasculature in contact with the proposed microenvironment cross-section is approximately 30 µ m. This latter distance is denoted as the cylindrical model's core, described as δ in Definition 3. Thus the radius of the cylindrical model is R = δ + ǫ , where ǫ is the aforementioned necrotic cuff. Refer to Fig. 1 to see a summary of spatial constraints. Using Proposition 1, subvolumes v i of cylindrical shells associated with compartments T i are calculated: The height of the cylindrical system h t , or level of tumor vascularization, is defined as a function with respect to time. The following linear vascularization function is assumed to simulate a slow-growing tumor: with k V denoting the rate of vascularization.
∀n ∈ N and basis case:  (3) carrying capacities to define constraints for their subpopulations. This subsection introduces cell motility and aggregate tumor outgrowth into more distant regions. For simplicity, cell metastatic migration is neglected in this process, and cells dying at each time step do not occupy physical space toward the capacity of the compartment. Additionally, some amount of cellular movement is allowed between the compartments to account for multi-directional cell proliferation and death. With a lack of data available for cell motility, the following probabilities (and resulting fractions) of compartment or subpopulation movement are used: where inward movement is slightly favored due to higher oxygenation. These probabilities were manually tuned via the tumor growth model performance. This multi-compartment system accounts for the self-interacting dynamics between subpopulations. Due to the addition of radial distance as a dimension, one must account for cellular growth that expands outward into a neighboring compartment, or oxygenation region, of the tumor. In the model, as cell populations reach their compartment's carrying capacities, they proliferate into more necrotic compartments under the assumption that the oncotic pressure decreases towards necrotic regions. Accordingly, the index notation i = 1 . . . n is employed to signify the ith tumor compartment T i and its corresponding proliferating population P i,t . Outgrowth into the i + 1th compartment from the neighboring ith compartment depends on the local spatial capacity of the ith compartment, which is denoted as K i,t . More precisely, the ith compartment's proliferating population must exceed that of its own capacity such that P i,t > K i,t . If this condition holds true, excess cells created in the current timestep t are henceforth contributing to the proliferating population of the i + 1th compartment, P i+1,t ( Table 2). This outgrowth (i.e., compartment overflow) at day t, O i+1,t , is defined as a recursive difference equation: This outgrowth term is added to P i,t+1 , compartment T i 's proliferative population at time t + 1 to yield a recursive definition of proliferative tumor population per compartment T i . With this spatiotemporal dynamism, each compartment is also governed by its corresponding local net growth rate k G,i : P(inward movement) := 0.10 P(outward movement) := 0.05 P(outward movement | 10-fold difference) := 0.10 protein shedding is represented as a relationship directly proportional to the number of cells in the system. Accordingly, the ith compartmental outflux at time t, u i (t) , is defined as a distance-based relationship between the number of tumor cells and the corresponding compartment's population and shedding values: where D i,t /�t is the dead tumor cell population for tumor compartment T i per day t, EC is the cellular localization of the protein ( EC = 1 denotes the extracellular domain and EC = 0 denotes the non-extracellular domain), is the protein-specific contribution per cell during necrosis, is the protein-specific shedding rate per cell per day, and w i is the distance-based weight for the ith tumor compartment, T i ( Table 1). Note that u i,t is the ith compartment's protein influx into the plasma at discrete day t. For the sake of simplicity, w i is approximated as a 1-dimensional diffusivity relationship with distance from vasculature: Finally, extracellular protein shedding is accounted for by active shedding from proliferating cells, while nonextracellular protein shedding is directly proportional to the instantaneous necrotic population at time t.
Heterogeneous plasma protein influx. Extracellular and non-extracellular protein mass in plasma, q(t) or q t , is modeled with the following ordinary difference and differential equations: where k E is the elimination rate of protein from the plasma, which encompasses protein degradation, excretion, clearance, etc. Equation (14) is approximated numerically as a difference equation in Eq. (15) using the Euler method. Parameter q 0 is the basal plasma protein mass, and relatedly, u H is the daily rate constant of influx of protein mass from healthy cells. It should be noted that when multiplied by t , the rate u H becomes an influx value at each time step. Furthermore, in experimental settings with tumor xenografts in mouse models, u H is assumed to be zero-valued since the generated human proteins of interest are non-endogenous. estimation of model parameter values. In order to build mappings between distance to vasculature r and net growth rate k G,i , data from both oxygen diffusion in tissue and in vitro tumor cell birth and death rates under varying oxygen conditions are used 20 . A compartment T i 's tumor growth parameters k G,i , k B,i , and k D,i is assumed linear with respect to oxygen concentration and estimated with the following equations: where C 0 C(0) is the volumetric concentration of oxygen at vasculature, r 1/2 is the distance half-life of oxygen, m is the slope describing the change in rates over the change in oxygenation, and b is the minimum measured rate value. Equation (16) assumes an exponential decay relationship 20 mapping radial distance from a proximal vasculature point r, to the volumetric concentration of oxygen C, thereby creating a concentration gradient over the domain r ∈ [0, ǫ] . Equation (17) assumes a linear relationship that maps the volumetric concentration of oxygen C to the net growth, birth, or death rate. Specifically, a linear interpolation is employed on birth and death rate data collected from cell cultures exposed to 1% and 20% oxygen to estimate n sets of birth, death, and net growth rates for the compartments midpoints s i as seen in Eq. (8) with the level set formulation. The resulting rates can be found in Table 1. A visualization of the birth, death, and net growth rates can be found in Fig. 2. Software implementation. The model was simulated using Python. Due to the coupled nature of tumor growth and protein shedding, these two phenomena were simulated in series, each with its own sets of parameters. The parameter set corresponding to tumor growth relies on the cell or tissue type, while the parameter set corresponding to protein shedding relies on the proteins of interest.

Source code availability
The source code for the model is available at https ://githu b.com/gmach iraju /Multi compa rtmen tTumo rs.
Received: 24 February 2020; Accepted: 10 September 2020 www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.