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.

Introduction

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-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-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 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 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 extra- and non-extracellular proteins from viable and necrosing cells as a function of heterogeneous vascularization.

Figure 1

Biological parameters for vascularized tumor growth and compartment model formulation. (A) 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 () 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 . 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 all generate cells with their specified birth, death, and net growth rates () 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 to account for heterogenous diffusion based on distance.

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 () and death rates ()[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.

Figure 2

Tumor compartment birth and death rates. Extrapolated rates are derived from cells exposed to varying levels of hypoxic stress[20]. The net growth rates () are equal to the birth rates () minus the death rates (). A nonlinear relationship can be observed as a function of compartment number (related to distance from proximal vasculature).

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 [24] has been measured in mammary adenocarcinoma[25,26]. This value can vary based on the tumor tissue type.

Necrotic cuff: The necrotic cuff () defines a region of viable cells and the space in which a gradient of cell proliferation exists. In situ, cells generally fail to grow beyond the diffusion distance of oxygen, approximately 100 m[21,23,27] from proximal vasculature. This metric generally holds for all tumor tissues.

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 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 protein-specific 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 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 () and the vascularization rate () were simulated. Parameter helps to determine the range of values for the compartmental net growth rates , 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 compartments, total cell population of greater than cells at 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.

Figure 3

Tumor growth dynamics over time and uniformly-selected simulations over a valid parameter space. (A) Selection of biologically-realistic tumors. This surface utilizes both of the parameters (vascularization rate) and (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 user-defined maximum iteration for simulation, . Each tile on the surface represents a simulation with the corresponding parameters on the x- and y-axes. If the value of is large and is too small, cells primarily stay in compartments near vasculature due to large carrying capacities and little outgrowth. If the value of is too small and 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 compartments (indicated by the color bar) and grow to a population of at least cells (represented by translucent plane) at the set maximum iteration of years. Tumors used for downstream simulation were both above the population plane and with 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 (B), (C), and (D). The trajectory that exceeded the arbitrary cutoffs while maintaining the smallest popuation size was used in all subsequent experiments and figures. (B) Valid simulations of tumor proliferation. Simulations are similar to that of Hori et al’s tuned model, but face a slowdown (i.e., greater cell death) in later stages due to unmet demands in vascularization. (C) Valid simulations of cell deaths per day. Daily cell death is largely caused by hypoxia-induced necrosis and begins nearly after years when the tumor reaches approximately cells. (D) Valid simulations of necrotic fraction represent the ratio of cumulative cell death to the total number of cells accounted for (cumulative cell death and current proliferative population). The necrotic fraction spikes starting at years when an uptick in cell death occurs at the same time, as seen in (C). This is an emergent property of the currently used model parameters, as is the case with the trajectories of (B) and (C).

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 and vascularization rate . Parameter 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 appears to enforce a ceiling function for the total tumor population after it outpaces vascularization and necrosis occurs at higher rates. For protein shedding, (protein half life in blood), , 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 and hyperparameters tuned to a particular tumor, with considerations such as cell and tissue type, while or and are parameters that can be estimated a priori for a shedding protein of interest.

Figure 4

Sensitivity analysis for tumor growth. Parameter values for (A) vascularization rate () and (B) the maximum volumetric concentration oxygenation at vasculature () interpolated were chosen at five uniformly-spaced values between their respective domain’s hypothesized minimum and maximum values, depending on the parameter of interest. The necrotic population trajectories are colored in black, while the proliferative population trajectories are colored in gray. For each plot, the parameter not undergoing analysis is fixed. Sensitivity of operated on a fixed value of =16 and varied values of . Sensitivity of operated on a fixed value of =0.1005 and varied values of . The summary of parameter values used for sensitivity are specified in Table 1, where fixed values are in the Value(s) Simulated column. It should be noted that boosting the value shifts the range of to more positive values, resulting in faster tumor growth. The inflection in the necrotic population (B) shows the uptick in cell death after the vascularization rate is outpaced by the growing cell population, causing more cells to experience hypoxic conditions.

Figure 5

Sensitivity analysis for protein shedding. Parameter values for (A) , (B) , (C) per-cell influx ( or depending on the population), (D) half-life , and (E) the healthy cell shedding influx , were chosen at five linear- or log-ordered points between their domain’s hypothesized minimum and maximum values, depending on the parameter of interest. The non-extracellular (non-EC) shedding trajectories are colored in black, while the corresponding extracellular (EC) shedding trajectories are colored in gray. For each plot, the parameter not undergoing analysis is fixed to a hypothetical set (=6.4, =0.00045 or =0.00045 depending on population, and =456) acting as baselines for both marker types. The summary of parameter values used for sensitivity are specified in 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 controls the immediate uptick in the trajectory, while parameters and or control the slope of the trajectory after proliferative growth slowdown and the emergence of the necrotic population. Parameter 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 and alone in the specified parameter ranges does not visibly benefit non-EC shedding since cell death is the source of markers.

Table 1

Model parameters

Parameter	Description (units)	Value(s) simulated	Range(s) simulated	References
Δt	Approximation time step (day)	1	-	-
k_G,i	Net tumoral growth rate of compartment T_i (day^-1)	See Eq. (17)	Calculated from k_B,i and k_D,i	–
k_B,i	Tumoral birth rate of compartment T_i (day^-1)	See Eq. (17)	[7.8e−4, 8.2e−3] as original range	[18]
k_D,i	Tumoral death rate of compartment T_i (day^-1)	See Eq. (17)	[0, 1.6e−3] as original range	[18]
P_0	Initial number of proliferating tumor cells (cell)	1	–	–
EC	Cellular localization of protein	{0, 1}	–	–
w_i	Distance-based diffusion weight for compartment T_i	See Eq. (13)	–	–
Ψ	Protein-specific contribution per cell (U/cell)	4.5e-4^#$	[1.4e−2, 4.5e2]^$, [4.5e−4, 4.5e−1]^+	–
Φ	Protein-specific shedding rate per cell (U cell^-1 day^-1)	4.5e−4^#, 4.5e−6^$	-	[8]
t_1/2	Blood half-life of protein (day)	6.4	[0.64, 6400]^+	[8]
k_E	Elimination rate of protein from plasma (day^-1)	–	Calculated from ln(2)/t_1/2	[8]
C_0	Maximum volumetric concentration of oxygen in the tumor (%)	16^&	[12,20]^∗+&	[18]
u_H,t	Healthy cell basal shedding influx; assumed constant (U/day)	4.56e2, 4.56e3^$	[4.56e1, 4.56e5]^+	[8]
k_V	Vascularization rate (day^-1)	0.101^&	[1e-3, 2e-1]^∗+&	–
p	Percent of total tumor volume occupied by tumor cells	0.2	–	[8]
ρ	Tumor cell density (cell/mm^3)	1e6	–	[8]
σ	Partitioning resolution, or width, of compartments (μm)	10	–	–
ϵ	Necrotic cuff (μm)	100	–	[15,19,20]
d	Mean vessel diameter (MVD) (μm)	60	–	[19]
r_1/2	Radial distance half-life of oxygen	0.018	–	[18]

U denotes arbitrary units. It is assumed that is multiple orders of magnitude larger than . Original ranges signify those found in the literature before any interpolation. Key: * Values/ranges simulated for valid tumor growth grid search (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).

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 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 and cells respectively, indicating that the simulated tumor reaches a steady-state in marker shedding at later stages.

Figure 6

Compartmental contributions toward tumor growth. Stack plot of the per-compartment (A) net proliferative and (B) necrotic populations over time. The tumor growth run used was programmatically selected for having both moderate tumor growth and necrosis, with the process of model selection seen in Fig. 1. The original exponential trajectory seen from 0–6 years slows down at approximately 6 years into the simulation due to vascularization failing to keep up with the growing metabolic demands of the tumor. This inflection is an artifact of the model’s assumed linear vascularization rate, i.e., the coupled mono-exponential growth models are limited by a linear ceiling. The majority of cells live and die in oxygenated compartments. However, there is a relatively larger contribution to cell death (B) from more distant compartments due to their very high death rates.

Figure 7

Compartmental contributions toward protein shedding. Stack plot of the per-compartment outflux (or vasculature influx) of (A) extracellular (EC) and (B) non-extracellular (non-EC) protein mass over time. The slowdown occurs due to the compartmental populations experiencing a limit in carrying capacity (based on vascularization) after the proliferative population reaches approximately cells. Shedding was simulated on the same tumor growth run used for Fig. 6, which was programmatically selected for having both moderate growth and necrosis. The parameter values of the simulated proteins were kept entirely identical, where specifically healthy cell influx (), and half-life () are the same. Furthermore, parameter values of were also assumed in this setting for ease of comparison. Once again, the relatively larger ratio of distant cell contribution to non-EC shedding (B) due to higher rates of cell death is visible here. The overall trajectory is near-linear (on log-scale) for both protein localizations.

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 , (healthy cell influx rate of protein mass), and . The colored regions show the dynamic range of trajectories for each value of , 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 ). 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.

Figure 8

Parameter scanning for non-extracellular proteins. Scans over protein parameters identifies cases when non-extracellular proteins outperform extracellular ones. Hypothetical extracellular (EC) and non-extracellular (non-EC) proteins are simulated and denoted by the lines embellished with diamonds and squares, respectively. Combinatorial scans over the parameters , , and are calculated. As discussed further in Fig. 5, parameter (set to values 4.56e+00, 4.56e+02, 4.56e+04) appears to control the immediate uptick in the trajectory, which can be interpreted as the initial levels of protein mass in circulation. Parameters (set to values 0.64, 20.24, 640) and (set to values 1.423e−02, 2.531e+00, 4.500e+02) appear to control the slope of the trajectory, which is especially visible after tumor growth slowdown and emergent necrotic population uptick. The minimum and maximum trajectories are taken for each value of scanned, resulting in shaded regions that define the operating dynamic range of shedding for the scanned parameter space of the other two parameter values. Specifically, the lavender-tinted, orange-tinted, and yellow-tinted regions are the dynamic ranges operating on varying the aforementioned and values along with set =4.56e+00, =4.56e+02, =4.56e+04, respectively. Due to a monotonic increase in protein mass with respect to and , given a set value of , the lower and upper bounds of dynamic ranges are always achieved with the smallest and largest values of (0.64 and 640) and (1.423e−02 and 0.423e−02), respectively. The crossover of dynamic ranges (tinted regions) indicates parameter operating regions where the initial boost in protein mass from is made up by the other parameters. Results are presented (A) with respect to time and (B) tumor population. The solid black line represents the Hori et al tuned shedding (using its underlying population growth model), while the line marked EC represents the shedding by the model proposed in this study, with the same parameter values as the Hori model. The end-behavior discrepancy between the two trajectories is most likely due to the differing tumor growth equations and associated assumptions.

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 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 population trajectories along with the corresponding shedding kinetics for a number of proteins of interest. 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 , 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 corroborates the need for density-dependent vascularization and faster tumor growth at later stages, as shown by the discrepancy of the presented model and the tuned model by Hori et al. This choice of 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, 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 (, , , and ) 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.

Methods

Proliferative and necrosing cell populations

To describe proliferation of a tumor cell population in vivo, a simple mono-exponential growth model is employed:where, P(t) is proliferative population at time and is net tumoral growth rate (Table 1). For model flexibility and cellular dynamics at each time step, difference equations are used at day intervals to form the Malthusian expression, . Through this use of a sufficiently small , Eq. (1)’s differential equation is converted to Eq. (2) via a forward Euler finite difference approximation:where is the net proliferative population at discrete time point and day is the time step taken for this approximation. The initial condition is set to to signify the number of tumor cells that initiate tumorigenesis. In order to approximate cell death at time t, the net tumoral growth rate is decomposed into its constituent birth rate and death rate :Using the above definition of net growth rate, the terms in Eq. (2) are expanded to define and , the amount of cell birth and death at time t, respectively:where 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, 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 , each with its own growth rate 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 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 of a real-valued function f of m variables is defined as a set of the form:where .

In order to parameterize the system, points in are mapped to cylindrical coordinates, giving rise to the spatial dimension of radial distance to proximal vasculature, r. Accordingly, the net tumoral growth rate is defined as a function of r:where the coordinate (x, y, z) is a point in and 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 tumor cell compartment () is assigned a level set that defines its corresponding volume and net growth rate . 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 is referred to as the local net growth rate:where is the radial midpoint of cylindrical compartment . The following mapping for is used to calculate all compartmental rate constants and shedding weights: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 of a cylinder is defined as the sub-cylindrical region starting at the origin with radius .

Definition 4

The concentric shell (corresponding to compartment ) of a cylinder has a volume of:

Definition 5

A cylinder with radius R and core is said to contain n equidistantly concentric shells if its radius is uniformly partitioned by . Namely, is the thickness of each cylindrical shell (i.e., associated with compartments ), where the concentric shell has radius:Furthermore, the shell subvolumes of s.t. 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 :

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 , where is the aforementioned necrotic cuff. Refer to Fig. 1 to see a summary of spatial constraints. Using Proposition 1, subvolumes of cylindrical shells associated with compartments are calculated:The height of the cylindrical system , 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 denoting the rate of vascularization.

Compartmental carrying capacities

The compartmental volumes help determine a compartment ’s associated physical capacity, . This can also be thought of as a numerical cell limit before cells begin outgrowth into the next compartment, i.e., being advected towards lower-pressure necrotic regions away from the vascular source by the oncotic pressure exerted by better oxygenated, proliferative cells. This relationship can be calculated by the following expression, using parameters measured from solid tumor cellular densities:where is the approximate tumor cell density (cells/mm) of a solid tumor, and p is the percent of total tumor volume occupied by tumor cells[8] (Table 1). Differing compartmental volumes account for differing carrying capacities. As briefly mentioned earlier with the introduction of a vascularization function , the volume and subsequent carrying capacity of a compartment changes over time. In effect, extends compartmental carrying capacities as blood vessels are recruited to increase tumor vascularization.

Compartment dynamics during tumor growth

The previous subsections fully defined tumor cell compartments with their own (1) effective rates of growth and shedding, (2) volumes, and (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 is employed to signify the tumor compartment and its corresponding proliferating population . Outgrowth into the compartment from the neighboring compartment depends on the local spatial capacity of the compartment, which is denoted as . More precisely, the compartment’s proliferating population must exceed that of its own capacity such that . If this condition holds true, excess cells created in the current time-step t are henceforth contributing to the proliferating population of the compartment, (Table 2). This outgrowth (i.e., compartment overflow) at day t, , is defined as a recursive difference equation:This outgrowth term is added to , compartment ’s proliferative population at time to yield a recursive definition of proliferative tumor population per compartment . With this spatiotemporal dynamism, each compartment is also governed by its corresponding local net growth rate :

Table 2

Model variables.

Variable	Description
P_i(t)	Proliferative population of compartment T_iat continuous time t (cell)
P_i,t	Proliferative population of compartment T_iat discrete time t (cell)
D_i,t	Dead cell population of compartment T_i at discrete time t (cell)
O_i,t	Cell overflow from compartment T_i to T_i+1 (cell)
h_t	Height of compartments (and cylindrical model) at discrete time t
K_i,t	Carrying capacity of compartment T_iat discrete time t (cells)
u_i(t)	Protein outflux from compartment T_i at continuous time t (U/day)
u_i,t	Protein outflux from compartment T_i at discrete time t (U/day)
q(t)	Protein mass in plasma at continuous time t (U)
q_t	Protein mass in plasma at discrete time t (U)

Heterogeneous intratumoral protein shedding

Consistent with previously established approaches[8], protein shedding is represented as a relationship directly proportional to the number of cells in the system. Accordingly, the compartmental outflux at time t, , is defined as a distance-based relationship between the number of tumor cells and the corresponding compartment’s population and shedding values:where is the dead tumor cell population for tumor compartment per day t, EC is the cellular localization of the protein ( denotes the extracellular domain and 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 is the distance-based weight for the tumor compartment, (Table 1). Note that is the compartment’s protein influx into the plasma at discrete day t. For the sake of simplicity, 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 non-extracellular 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 , is modeled with the following ordinary difference and differential equations:where 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 is the basal plasma protein mass, and relatedly, is the daily rate constant of influx of protein mass from healthy cells. It should be noted that when multiplied by , the rate becomes an influx value at each time step. Furthermore, in experimental settings with tumor xenografts in mouse models, 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 , 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 ’s tumor growth parameters , , and is assumed linear with respect to oxygen concentration and estimated with the following equations:where is the volumetric concentration of oxygen at vasculature, 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 . 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 and oxygen to estimate n sets of birth, death, and net growth rates for the compartments midpoints 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.