Using Fractal Geometry and Universal Growth Curves as Diagnostics for Comparing Tumor Vasculature and Metabolic Rate with Healthy Tissue and for Predicting Responses to Drug Therapies

Healthy vasculature exhibits a hierarchical branching structure in which, on average, vessel radius and length change systematically with branching order. In contrast, tumor vasculature exhibits less hierarchy and more variability in its branching patterns. Although differences in vasculature have been highlighted in the literature, there has been very little quantification of these differences. Fractal analysis is a natural tool for comparing tumor and healthy vasculature, especially because it has already been used extensively to model healthy tissue. In this paper, we provide a fractal analysis of existing vascular data, and we present a new mathematical framework for predicting tumor growth trajectories by coupling: (1) the fractal geometric properties of tumor vascular networks, (2) metabolic properties of tumor cells and host vascular systems, and (3) spatial gradients in resources and metabolic states within the tumor. First, we provide a new analysis for how the mean and variation of scaling exponents for ratios of vessel radii and lengths in tumors differ from healthy tissue. Next, we use these characteristic exponents to predict metabolic rates for tumors. Finally, by combining this analysis with general growth equations based on energetics, we derive universal growth curves that enable us to compare tumor and ontogenetic growth. We also extend these growth equations to include necrotic, quiescent, and proliferative cell states and to predict novel growth dynamics that arise when tumors are treated with drugs. Taken together, this mathematical framework will help to anticipate and understand growth trajectories across tumor types and drug treatments. 1. Introduction. As tumors grow, they must develop a vascular system to supply an ever increasing need for oxygen and nutrients. Through the process of angio-genesis, vessels permeate the tumor and deliver these oxygen and nutrients through blood[1, 16, 28, 58, 74]. These resources are used to fuel the maintenance and growth of tumor cells. The tumor vascular system appears to be relatively autonomous when removed from the body. Indeed, measurements of vascular structure and blood flow are often taken for this isolated system by perfusing the tumor with blood and interpreting the findings in isolation of the host[55, 58, 81]. Very little quantification of tumor vascular properties exists have been proposed to capture the structure, variation, or mechanistic principles that characterize tumor vasculature, calling for a larger role of mathematical modeling in studying solid tumors [22, 10, 9]. Considering tumor vascular networks in isolation greatly simplifies the modeling approach and allows existing …

Abstract.Healthy vasculature exhibits a hierarchical branching structure in which, on average, vessel radius and length change systematically with branching order.In contrast, tumor vasculature exhibits less hierarchy and more variability in its branching patterns.Although differences in vasculature have been highlighted in the literature, there has been very little quantification of these differences.Fractal analysis is a natural tool for comparing tumor and healthy vasculature, especially because it has already been used extensively to model healthy tissue.In this paper, we provide a fractal analysis of existing vascular data, and we present a new mathematical framework for predicting tumor growth trajectories by coupling: (1) the fractal geometric properties of tumor vascular networks, (2) metabolic properties of tumor cells and host vascular systems, and (3) spatial gradients in resources and metabolic states within the tumor.First, we provide a new analysis for how the mean and variation of scaling exponents for ratios of vessel radii and lengths in tumors differ from healthy tissue.Next, we use these characteristic exponents to predict metabolic rates for tumors.Finally, by combining this analysis with general growth equations based on energetics, we derive universal growth curves that enable us to compare tumor and ontogenetic growth.We also extend these growth equations to include necrotic, quiescent, and proliferative cell states and to predict novel growth dynamics that arise when tumors are treated with drugs.Taken together, this mathematical framework will help to anticipate and understand growth trajectories across tumor types and drug treatments.1. Introduction.As tumors grow, they must develop a vascular system to supply an ever increasing need for oxygen and nutrients.Through the process of angiogenesis, vessels permeate the tumor and deliver these oxygen and nutrients through blood [1,16,28,58,74].These resources are used to fuel the maintenance and growth of tumor cells.The tumor vascular system appears to be relatively autonomous when removed from the body.Indeed, measurements of vascular structure and blood flow are often taken for this isolated system by perfusing the tumor with blood and interpreting the findings in isolation of the host [55,58,81].Very little quantification of tumor vascular properties exists (but see [23,55,56,58,75]), and until recently, few models [1,2,11,89,31] have been proposed to capture the structure, variation, or mechanistic principles that characterize tumor vasculature, calling for a larger role of mathematical modeling in studying solid tumors [22,10,9].Considering tumor vascular networks in isolation greatly simplifies the modeling approach and allows existing frameworks for healthy tissue [101,102,103,20,97] to be straightforwardly applied to tumor vascular systems.This application allows direct comparison with vascular structure in healthy tissue and assessment of whether similar principles are at work in the development of both systems (see also [72]).
Leading models for the structure of the healthy vascular system, such as the West, Brown, Enquist (WBE) model [98,76], rely on just a few key principles.First, it is posited that the network is hierarchical and branching.This assumption allows the network to be described by scaling ratios for vessel numbers, radii, and lengths across levels k in the network.Second, energy minimization is proposed to constrain the scaling ratio for vessel radii.This constraint results in area-preserving branching (for the total cross section across levels) for large vessels in which pulsatile flow and wave reflections dominate and results in area-increasing branching (i.e., Murray's Law [67]) in small vessels for which dissipation and viscous forces dominate.Third, the network is assumed to be space filling [61] to ensure capillaries are close enough to each cell to supply oxygen by diffusion, thus placing a constraint on the scaling ratio for vessel lengths.Fourth, capillaries-the terminal unit of the network-are recognized as having properties that are invariant with body size or tissue type [77].As a result, the capillary properties set the scale for the entire network.This last assumption is the one most easily adapted to tumor vascular structure.As explained below, all of the other assumptions may be less well grounded in the principles that drive tumor growth and development.Nevertheless, the scaling rules for vessel branching in healthy tissue are hypothesized to be a result of evolution and have been shown experimentally to be encapsulated in a genetic program for vessel growth [26,65].Tumor growth likely co-opts these existing genetic programs but may introduce substantial noise into their execution and expression.Through theory and data, we can begin to assess the degree to which this holds.
Despite the simplifying assumptions above for an autonomous network, within the body the only source of blood is from the heart.Therefore, the tumor vascular system latches onto one or more vessels in the healthy cardiovascular system and diverts flow from healthy tissue to the tumor.Necessarily then, there is an intimate integration between the host and tumor vascular networks.We incorporate some of these integration effects by using existing results from models for the vascular density and blood flow rate in healthy vasculature.This allows us to determine the likely number and flow rates of vessels that serve as sources for the blood flow in the tumor.
The growth of solid tumors beyond small sizes depends on oxygen and nutrients supplied initially by diffusion from nearby capillaries from healthy tissue in the host [46,47].Beyond small sizes, however, growth can only occur through the process of angiogenesis and vasculogenesis within the tumor.Metastasis allows cell from a primary tumor to bud off and establish new tumors in other locations that then repeat this process of oxygen supply via diffusion and eventually angiogenesis [59,69].Again building on previous results for universal ontogenetic growth curves across species [98,24,25,14], we reinterpret parameters in terms of tumor characteristics to derive growth curves for tumors [4,15,48,68,83,91,92,93] that can be adapted to specific tumor types.
These models and findings are of critical importance because tumor vasculature has been targeted for treatment in a variety of cancers and models [49,50,51,52,53,54,88,73]. Early on, Folkmann [16] and others proposed the idea of stopping formation of the vascular system within the tumor and starving the tumor.More recently, researchers have been studying a different approach that involves remodeling tumor vasculature to more-closely resemble healthy vasculature and thus deliver blood more uniformly throughout the tissue, enabling better drug delivery [1,2,43,42].
However, even with improved drug delivery most treatments for cancer face an additional obstacle: most cancer drugs target only proliferating, (i.e., actively dividing) cells.Within tumors, however, many cells become quiescent such that they are still alive and metabolically active but are no longer dividing [12,18,62,63,7,95,30,21].After the drug concentration in a tumor drops, these quiescent cells can transition to a proliferative state, thus evading treatment and allowing the tumor to grow back.We outline a model that identifies spatial gradients and thresholds within the tumor for delivery of energy and resources.For large tumors, the core becomes oxygen-and nutrient-poor and thus necrotic [3].For annuli of tissue at intermediate radii, cells will receive enough oxygen and resources for basic maintenance but not for replication (i.e., quiescent state), and at larger radii enough power is available for both maintenance and growth (i.e., proliferative state).
This work is the first to use a single framework to couple the mathematics of fractal networks with the temporal and spatial growth of tumors.Our paper is organized as follows.We first derive the equations for the scaling of vascular networks and use geometric properties of these networks to derive metabolic rate.Next, we combine these results for metabolic rate with universal growth equations to predict growth curves for tumors.Third, we expand the growth equations to consider the spatial distribution of different metabolic states of cells, corresponding to different annuli within the tumor.Finally, we compare to known empirical results and consider the implications of our results for the mechanisms that control tumor vascular structure and growth.

2.
Structural and scaling properties of vascular networks.In this section, we review the basic framework for healthy vascular systems, discuss how these principles each apply to tumor vasculature, and discuss how the framework and principles may need to be adapted for future models.
2.1.Notation for branching, hierarchical network.For healthy vascular systems, it is assumed that the average properties of the network can be largely captured in terms of a symmetric, branching, hierarchical network [97,76].Specifically, the network begins with a heart or feeding vessel through which blood flows through successive branching junctions until it reaches the capillaries that deliver oxygen to the cells.At each branching junction, a parent vessel branches into n identical daughter vessels.The average branching ratio is applied at every branching junction in the model version of the network such that where N k is the number of vessels at level k and n does not vary within or across levels k.This branching ratio transcends local and global scales in the network and thus becomes the key parameter that relates all of the other variables and characteristics of the network.
Figure 1.Plot of the logarithm of the number of vessels.N k versus the branching order, k, as determined by the Horton-Strahler method [34,87].Data are from [37,99,100,45,35] from 8 healthy sections (e.g., left anterior descending artery, pulmonary) of vascular systems in pig, cat, and human.In all cases, the correlations are excellent and corresponding branching ratios are near n = 3. See Fig. 6 for estimates of scaling exponents based on best fit regression lines as well as p-values, 95% confidence intervals, and r 2 values.
Figure 1 shows extremely systematic branching of vasculature from healthy tissue in mammals.The average branching ratio for these data is closest to the integer value of n = 3. Visual inspection of vessel casts and images suggest bifurcating (n = 2) branching, but the discrepancy is likely due to the choice of Horton-Strahler ordering for measuring these data.The Horton-Strahler method labels all terminal vessels as the same level, and the index of upstream branches is determined by the sum of the two downstream branches that connect with it [34,87].Consequently, a branching ratio of n = 3 means that there are three times as many vessels in one level compared to the level that feeds into it because of this choice of labeling with the Horton-Strahler method.For any specific branching junction, however, it does not mean that there are three daughter vessels.Because all empirical data collected thus far have been measured using the Horton-Strahler method, we must follow that notation here to be consistent for all calculations, and thus, we will assume a branching ratio of n = 3 throughout the rest of the paper.We are in the process of automating measurements for the vascular system that do not use the Horton-Strahler method, but instead label the level of a vessel by how many branching junctions it is away from the heart [97,76].
Within healthy tissue, there are violations of hierarchical, symmetric branching [101,103].For example, in some parts of the vascular system, such as the coronary artery, daughter vessels are often not identical at a branching junction, resulting in asymmetric branching.Moreover, in other parts of the body, such as the circle of Willis in the brain, the network is known to contain loops that break the strict assumption of hierarchy and lead to more grid-like geometry.For healthy vasculature, however, the symmetric, hierarchical network seems to do a reasonable job of capturing the basic features of the network.
For vascular systems within tumors, images suggest the network contains many more loops and is much more grid-like than for healthy vasculature.Consequently, including new metrics for network structure may prove important as continued progress is made in modeling tumor vascular networks.In particular, some very recent models for reticulate venation networks in leaves with high levels of loopiness may help give insight into the best way to model tumor vasculature [70,8,36].For now, we proceed with models for healthy vascular systems in order to facilitate comparisons [97,76].

2.2.
Network branching that spans body to deliver blood to all cells: Constraint on vessel lengths.In order for a capillary to be near enough to each cell that oxygen can be delivered by diffusion, each vessel at a given level must itself supply a volume of smaller vessels (ultimately capillaries supply a volume of cells) in such a way that the total service volume is conserved from level to level.This service volume is approximated by a sphere (or ellipse, rectangle, or other Euclidean shape) with radius proportional to vessel length l k .The sum over these service volumes within a level k must equal the whole volume, which is the body that it supplies.Applying this rule across successive branching levels and service volumes gives a constraint on vessel lengths.
Because the right side of this equation does not depend on level k, the relative lengths of vessels from one level to the next is the same, independent of location in the network or absolute size of the vessels.Thus, the space-filling constraint exactly predicts a specific form (−1/3 scaling) of a self-similar fractal for the cardiovascular system.Specific predictions of fractal forms and specific scaling exponents that are based on physical or biological principles are rare and help provide a firmer footing for the prevalence of fractals in the natural world.Because of necrosis, the vascular system for tumors does not span the entire volume of the tumor.Instead, within the tumor there are large gaps with no vasculature that result in necrotic patches [3].We assume, however, that the tumor vascular system spans its entire viable volume.
2.3.Minimization of power to deliver blood: Constraint on vessel radii.Power is predominantly lost in pumping blood through the vascular system for two reasons.First, blood is lost due to reflections of pressure waves at branching junctions.These reflections can be completely eliminated using impedance matching with impedances derived from the Navier-Stokes equations in the large vessel limit for flow through elastic tubes [97,76].This derivation results in area-preserving branching for vessels with large radii such that As for vessel lengths, the right side of this equation is independent of level k, once again implying a very specific form of a self-similar fractal for the vascular system.Second, power can be lost due to dissipation and viscous forces between the blood and the vessel walls.This source of power loss is most important for small vessels with large surface-area-to-volume ratios.This type of power loss is non-zero for any vessel with a finite radius.Thus, minimization of this power loss must be done subject to constraints using a Lagrange-multiplier calculation [97,76].The constraints imposed are that the network volume must fit within the body and the network must span the space to feed all cells (section 2.2).The result is areaincreasing branching Fractal behavior is hence predicted once again, and for vessel radii, there are two regions with approximate self-similar fractal scaling: large vessels with area-preserving branching and small vessels with systematic area-increasing branching.
For tumor vasculature, minimizing power loss may still be advantageous because it conserves energy available for maintenance or growth.If the scaling principles for healthy tissue are analogous, then most tumors would be expected to exhibit area-increasing branching because their vessels will primarily be in the small vessel regime.On the other hand, tumors might be more interested in maximizing available power and thus maximizing growth to evade the immune system and drugs.Thus, some other optimization or joint optimization might lead us to expect different values for β or a deeper understanding of its observed variations.
2.4.Generalization of scaling ratios and exponents to arbitrary tumor vascular networks.The scaling exponents predicted and measured for vascular networks in healthy tissue serve as a useful baseline and point of departure for tumor vascular networks.Nevertheless, as mentioned above, the design principles controlling the structure and dynamics for healthy vasculature may not apply to tumor vasculature.Consequently, we maintain the assumption of a vascular network that on average can be represented as a self-similar fractal, but we relax the strict assumptions and predictions from existing models by allowing the specific value of the scaling exponents to be unconstrained and measured directly from empirical data.Such an approach has previously been done for plants [71,78,31].In particular, we let the scaling ratios be expressed as and represent the average scaling properties of the network.Because of the error distribution for vessel radii and lengths, it is most appropriate to measure these exponents in logarithmic space [79,39].Moreover, because vessel radii and lengths are often measured across levels k, it is easiest to express these scaling ratios relative to some fixed vessel (terminal or feeding).By repeated application of the scaling ratio, we derive and where l 0 and r 0 are the length and radius of the feeding vessel that supplies blood for the whole network.In logarithmic space, these expressions become Thus, by plotting ln(l k ) versus k and ln(r k ) versus k and knowing the value of the branching ratio n, we can measure the scaling exponents b and a from the slopes.Figures 2-5 plot available data for the scaling ratios of lengths and diameters (or radii) for both healthy and tumor vasculature.Tumor vasculature differs from healthy in at least two notable ways.First, the slopes (scaling exponents) are shallower for both diameter and length scaling.These results suggest vessels in successive levels are closer in size and that the terminal vessels are larger in tumors.Second, the tumor data exhibits more curvature around the line, worse fits, and quite large variation.These noisy patterns could be a result of the higher mutation rate and chromosomal instability in tumor cells that introduce noise into the genetic machinery for vessel development [26,65].Together, these two results indicate a breakdown in hierarchy for tumor vasculature and much larger role for variation in scaling and size.Hence, new types of models may need to be developed that can better incorporate these factors and describe these novel aspects.
Values for the scaling exponents are calculated in Fig. 6 based on the measurements from Figures 2-5 and a branching ratio of n = 3 based on Figure 1.Once again, the scaling exponents are much shallower and exhibit larger variation than is observed for healthy tissue.This conclusion holds when comparing either means or individual tissues to individual tumors.
The values in Fig. 6 can be used as input for the equations derived in the section below to predict how blood flow rate and metabolic rate scale with tumor size for various tumor types.
3. Blood flow and metabolic rate predicted from geometry of vascular structure.With knowledge of the scaling ratios for vessel radii and lengths along with the size of the capillaries to set the overall scale, we can now calculate many structural properties (e.g., surface areas and volumes for vessels or entire levels and subsets of the network) as well as flow rates and impedances for vessels, entire levels, or network subsets.As an example, we calculate the volume of the entire network and relate it to the number of terminal units for the case when a and b are fixed across the entire network.Assuming all vessels are cylindrical in shape and using Data are from [37,99,100,45,35,56,55] from 7 healthy sections (e.g., left anterior descending artery, pulmonary) of vascular systems in pig, cat, and human.For healthy tissue, the correlations are excellent and correspond to scaling exponents that on average are close to the value of b = 1/2 (see Fig. 6), somewhat higher than predicted.See Fig. 6 for estimates of scaling exponents based on best fit regression lines as well as p-values, 95% confidence intervals, and r 2 values.scaling equations similar to Eqs. (5)(6) to simplify our sum, the total volume of the network is given by summing the vessels across all levels.
where N N is the total number of capillaries and V N is the volume of a single capillary, which is invariant according to one of our primary assumptions.Substituting in Eqs.(7)(8) yields When N N is very large and 2a + b > 1, this equation can be approximated as Figure 3. Plot of the logarithm of the length of vessels, l k , versus the branching order, k, as determined by the Strahler method.Data are from tumors (mammary and colorectal carcinomas (A91 and A97)) in rat and human.For tumors, the correlations are much weaker and the scaling exponents are significantly shallower than for healthy tissue shown in Fig. 2. For both healthy tissue and tumors, a self-similar fractal network is a reasonable approximation for the vascular systems, but for tumors there is still some significant variation to explain.See Fig. 6 for estimates of scaling exponents based on best fit regression lines as well as p-values, 95% confidence intervals, and r 2 values.
Based on Eqs.(2)(3) this predicts that for large vessels and ) for small vessels or also for cases when 2a + b < 1.Data are from [37,99,100,45,35,56,55] from 7 healthy sections (e.g., left anterior descending artery, pulmonary) vascular systems in pig, cat, and human.For healthy tissue, the correlations are excellent and correspond to scaling exponents that on average are close to the predicted value of a = 1/2 for larger vessels (see Fig. 6).Also see Fig. 6 for estimates of scaling exponents based on best fit regression lines as well as p-values, 95% confidence intervals, and r 2 values.
Note that for the size ranges observed across species, finite-size corrections can significantly impact the exact form of these relationships [76,79].Over the range of sizes for tumors, which is much smaller than body sizes across species, these effects may be especially pronounced, but for the zeroth-order model presented in this paper, we ignore these finite-size corrections to first order.
Finally, if we assume that all capillaries are identical in terms of blood flow rate Q and metabolic rate, B, produced from the supplied oxygen, and further note that network volume scales directly with body volume and mass, M , [97,76] we derive predictions for these dynamical properties as well.We find for large vessels when 2a + b > 1 or for small vessels or also for cases when 2a + b < 1.In all cases, the values of the fractal scaling exponents, a and b, play the dominant role in determining the scaling of metabolic rate with body mass.Data are from tumors (mammary and colorectal carcinomas (A91 and A97)) in rat and human.For tumors, the correlations are weaker and the scaling exponents are significantly shallower than for healthy tissue shown in Fig. 4. For both healthy tissue and tumors, a self-similar fractal network is a reasonable approximation for the vascular systems, but for tumors there still remains significant variation to explain.See Fig. 6 for estimates of scaling exponents based on best fit regression lines as well as p-values, 95% confidence intervals, and r 2 values.
4. Universal growth equations.Understanding how tumor metabolic rate and blood supply change with tumor size as well as with properties of the host and specific tissue [96,31] will enable us to derive equations for the growth trajectory of the tumor [4,15,48,68,83,91,92,93].Elaborations of this first simple model will also permit us to include effects of multiple metabolic states (Sec.5) and eventually include effects of drug delivery and pressure (Sec.6).As a first step, we employ an equation that requires conservation of energy [97] and makes predictions for a universal ontogenetic growth curve that well matches empirical data.Guiot et al. [24,25] have previously studied the universal properties of tumor growth, but they did not interpret their fitted parameters in terms of the geometry of the tumor vascular system or host properties.We have explored universal functional forms for tumor growth [31] in a previous publication, but we did not compare to ontogenetic growth data or include effects of quiescent cells or treatment by drugs.
At the most basic level, oxygen and nutrients are used to power maintenance of existing cells and growth of new cells through replication and division.Thus, the total metabolic rate, B T , of the tumor can be expressed in equations as two terms, the first of which is power devoted to maintenance of existing cells and the second of which is devoted to the overall growth of the tumor through the creation of new cells where N v is the number of viable cells at time t, N c is the total number of cells, B m is the power devoted to maintaining each cell, and E c is the energy to create a cell.The total number of viable cells increases according to the net difference between the rate of cell birth and death: where N d is the number of cells that have died by time t.
At steady state (i.e., dN v /dt = 0), cells are still dying and being replaced.Following Herman et al., we first assume that each cell has the same chance of dying, independent of tumor size or location in tumor.Therefore, we can define the inverse lifetime of an average cell, Γ, by In real tumors, this inverse lifespan will almost certainly change with location in the tumor and also tumor size due to differences in access to resources (i.e., highly vascularized versus non-vascularized regions) and higher pressures near the center of the tumor [5,84].Combining Eqs. ( 19) and ( 20) we can re-express the growth equation, Eq. ( 18), purely in terms of viable cells: where we define B c ≡ B m −ΓE c to be the average cellular metabolic power required for maintaining a tumor at steady state.Because the viable tumor mass is , where N v is the number of viable cells and it follows that m c is the average mass of a tumor cell, we can further simplify our growth equation as Isolating the time derivative of the viable mass to one side of the equation yields Using our results from the previous section, we replace the tumor metabolic rate, B T , with its dependence on viable (metabolically active) mass via the network geometry and measured fractal scaling exponents, a and b, for 2a + b > 1 and for 2a + b < 1. B 0 is a normalization constant that captures the dependence on host properties by relating tumor metabolic rate to tumor size.This normalization contains the overall blood supply rate of the tumor via the feeding vessels as determined by the number and flow rate of nearby vessels in the host.This first-order differential equation now explicitly links properties of tumor cells (B c , E c , and m c ) with properties of the whole tumor (B T and m v ), the geometry of the vascular network (a and b), and those of the host (B 0 ).
For this baseline model, the equation can be solved exactly in either case.When 2a + b > 1, we solve to obtain where M v is the asymptotic viable mass of the tumor and D = m c B 0 /E c .When 2a + b ≤ 1, which is naively expected for early tumor growth driven by diffusion or the small vessel regime that minimizes energy loss to dissipation, the solution reduces to the especially easy form of exponential growth with α = (m c B 0 − B c )/E c determined by both cellular and host properties.Intriguingly, for all of the tumor data we analyzed, 2a + b ≤ 1, implying exponential growth.
In later stages of growth, tumors are expected to have large vessels that are potentially capable of pulsatile flow [97,76], and may be expected to alter the geometry such that 2a + b > 1, resulting in classic sigmoidal growth with the viable mass reaching a fixed asymptotic value given by and is attained when t ).As shown in Herman et al. [31], these predictions for growth trajectories are in good agreement with available data.This is impressive given the biological detail ignored by this model.
It is important to note that we consider each tumor in isolation and do not distinguish between primary and metastatic tumors or account for effects of metastasis at all [59,69].When tumors reach an asymptotic mass as predicted above, or even before this happens, metastasis may occur, thus enabling tumor cells to continue growing in the host by finding new locations and new local resources to draw upon.Indeed, this process of metastasis is likely the primary method by which tumor cells can circumvent an asymptotic size, continue growth, and eventually lead to more global failures for the host.
The general solution for the growth equations can be made dimensionless in order to predict a universal growth curve [98,24].Specifically, dimensionless mass is defined by where M is the asymptotic mass of the tumor.Dimensionless time is defined by where m 0 is the initial or birth mass and M is asymptotic mass.Using this equation, fits can be performed to determine the unknown parameters, and data from all growth trajectories should collapse to a single curve.
Figure 6 shows data for ontogenetic growth trajectories from birth to adult across diverse species [98].Data for tumor growth trajectories from initial cells to death of host are also plotted [31,92].Remarkably, these two types of growth (tumor and ontogeny) do approximately collapse to a single common curve.This is despite the fact that these two processes of growth are dramatically different in process, cell differentiation, cell regulation, and many other properties.Furthermore, this finding indicates the advantage of studying growth in terms of energetic costs and supply rate, providing a common currency for growth.Finally, we note that the tumor  29)- (30).Data in grey are for ontogenetic growth from 13 species of animals (see [98] for original data sources), ranging from guppy to cod to guinea pig, and data in color with square symbols are for tumor growth trajectories for C3H mammary carcinoma (dark blue), EMT6 mammary carcinoma (dark red), KHJJ mammary carcinoma (light green), NCTC (dark green), Flank (yellow), Primary fibroadenoma (red), Primary Osteosarcoma (light blue), and Walker Carcinoma (purple) tumor types (see [31,92] for original data references).This plot shows the extent to which all the growth data can be described by a single universal curve [98,24], and thus, the extent to which tumor growth is really analogous to ontogenetic growth form birth to adult.Fitted parameter values are taken from the original references.growth trajectories are dominated by the exponential phase of growth as compared to the ontogenetic curves, and that this phenomenon is consistent with the shallower scaling exponents observed for tumors and the condition that 2a + b < 1.The tumor growth trajectories do exhibit some slowing in growth at late time, but they never show the clear asymptote observed in ontogenetic growth.Perhaps tumor vasculature always enables exponential growth until death of the host, whereas ontogenetic growth and vasculature must eventually change and asymptote to enable survival and maintenance of existing cells.This suggests that tumor growth matches most closely with early ontogenetic growth, and that understanding early ontogeny may be most relevant to helping understand tumor growth.Conversely, further study of how ontogenetic growth changes from early to late ontogeny may help give insights into possible tumor treatments.Further theoretical and empirical studies are needed to illuminate these and other questions about the growth of tumors.

5.
Coupled growth equations for cells in proliferative, quiescent, and necrotic metabolic states.We now show how this model can be improved by incorporating cells of different metabolic states [12,18,62,63,7,95,30,21].These differences in metabolic state likely arise due to differences in resource supply and pressure varying throughout the tumor [5,84].We allow for the possibility of proliferative (i.e., actively dividing) cells, quiescent (i.e., metabolically sustaining but not dividing) cells, and necrotic (i.e., dead but not yet removed) cells.Each cell type requires an equation, and the equations are coupled by transition rates between each of the cell types.
where m is the mass, p is proliferative cells, q is quiescent cells, n is necrotic cells, f p→q is the frequency of conversion of p to q, Γ q is the rate of necrosis, λ is the rate of lysis for necrotic cells, S T is the surface area of the tumor, and τ M is the endothelial cell density for the host.In principle and often in practice, parameters such as f p→q and Γ can depend on tumor size, metabolic rate, and pressure.Only more explicit mechanistic models prescribe what those dependencies are.In these equations, we have made the approximation that the viable or metabolically active mass of the tumor scales linearly with the proliferative mass.This assumption is reasonable because the proliferative cells have much higher metabolic rates and are usually much larger in numbers, allowing them to dominate the dependence of metabolic rate on size within the tumor.During early stages of growth, nearly all of the tumor cells are proliferative, so this approximation will strongly hold for the early part of the growth trajectory.
Note that proliferative cells can only transition to become quiescent (f p→q ) and are not allowed to effectively skip the quiescent stage and become necrotic.Furthermore, necrotic cells cannot transition to any other metabolic state because there is no resurrection of dead cells.In general, these equations cannot be solved analytically, so numerical methods must be used.It is straightforward to obtain numerical solutions and growth trajectories using Mathematica.Figure 7 shows an example growth trajectory for these three types of cells within the tumor.
This model is based on a fundamental difference from growth in healthy tissue where all three metabolic states do not typically occur within the same cell type or tissue.
5.1.Thresholds for metabolic states at different layers within the tumor as driven by differential access to resources.The equations above are generally applicable to any coupled set of cell types that can transition from one to another.The biological and physical principles specific to tumors enter through the values of the parameters, many of which are similar to those given above and in Herman et al. [31].The new parameters here are the transition probabilities, f s, and the clearance rate λ.
A region of necrotic cells is often formed at the core of tumors due to lack of oxygen (hypoxia) [3].This hypoxia can occur because vessels become closed due to higher pressures and the vasculature delivers the most blood to the outer regions of the tumor.Similarly, a region of quiescent cells may form at an intermediate radius in the tumor, representing a transition region from the most actively dividing proliferative cells at the outer edge to the inner necrotic core [12,18,62,63,7,95,30,21].
Here, we show how to incorporate a simple radial dependence into the transition rates in the growth equation Eq. ( 33) by assuming available oxygen decreases as a function of distance from the edge to the center of the tumor.We place thresholds on the supply rate of metabolic energy necessary to sustain proliferative and quiescent cells, thus determining the transition rates as a combined function of radius and overall tumor size.The location of this threshold will continually change based on the growth of the tumor and the ability to access feeding vessels in the host.
The metabolic rate of the tumor is B T and is posited to be proportional to the tumor surface area, S T , because it is this area through which feeding vessels can be recruited.Thus, the total metabolic rate of an annulus of cells at radius r is B in , so where k(r) is the fraction of the total metabolic rate for all of the cells in an annulus at radius r with thickness equal to the diameter of a single cell, dr, and R is the radius of the whole tumor.Because B T does not depend on r, To perform further calculations, we must assume a shape for the tumor so that a functional form for its surface area, S T is specified.Here, we assume the overall tumor shape is approximately spherical.More generally, we are assuming a Euclidean geometry (as opposed to a fractal geometry like the tumor vasculature) for the overall shape of the tumor.There is no inherent conflict in assuming a fractal geometry for the scaling of vessel dimensions across branching levels in the vascular network and simultaneously assuming a Euclidean geometry for the shape of the whole tumor.We make this assumption of a Euclidean shape because to our knowledge there is no general evidence for a non-Euclidean or fractal surface for the whole tumor.If the overall tumor shape is indeed fractal, this can be straightforwardly incorporated into our modeling framework by using a different functional form for S T in all of the equations.
We assume that metabolic intake scales as a power law such that k(r) = cr α , with a normalization constant c and exponent α .Using the normalization condition to solve for c yields the expression c = α +1 R α +1 .Because we assume the tumor is spherical, the tumor surface area is 4πR 2 , so where K is the normalization constant that relates tumor surface to total blood and power supply to the tumor via feeding vessels.The size and number of feeding vessels are in turn determined by host properties, just as B 0 was in Sec.(4).The metabolic rate of each cell, B c (r), is calculated by dividing B in (r) by the number of cells in the annulus at radius r with thickness dr equal to the diameter of an average tumor cell.
• For α < 2, B c (r) would be inversely proportional to the radius.This would indicate that the innermost cells are the most active, which is the reverse of what has been observed experimentally.• For α = 2, B c (r) would not depend on r.
• For α > 2, B c (r) will be directly proportional to a power of the radius.
The transition point at α = 2 makes intuitive sense because it represents the point at which the change in resources is matched with the change in the number of cells per annulus across the surface area.
We set B c (r) equal to B c,q , the power below which proliferative cells will enter the quiescent state to determine the threshold for the largest radius at which proliferative cells can no longer exist.r outer,q = r inner,p = 3 4πKdr 2 (α + 1) B c,q Similarly, we set B c (r) equal to B c,n , the power below which cells will die and become necrotic, in order to find the largest radius at which quiescent cells become necrotic, or inversely the smallest radius at which quiescent cells exist.We find This calculation shows that, given α > 2, the threshold radii for quiescence and necrosis will increase faster than the radius, R, of the whole tumor.This effect will continue until the proliferative cells in the outermost annulus of the tumor become quiescent, i.e., the tumor can no longer grow due to environmental constraints.When the tumor is approximately spherical in shape, it also holds that where ρ T is the density of the tumor.This expression follows because the outer radius of the tumor, R, is the only variable that depends on time, and the sum of the necrotic and quiescent volumes must equal the volume contained within the radius of the outer layer of quiescent cells.Consequently, in this scenario, tumors must reach an asymptotic size, and exponential growth cannot continue forever, as was potentially possible for the growth equations that incorporated only proliferative and necrotic cells in Sec.(4).Moreover, the growth model presented in Sec. ( 4) and Herman et al. [31] predicts necrosis but not its spatial location.Indeed, for uniform vasculature, which is the tacit assumption in those growth models, the necrosis will most likely occur at the outer layer of the tumor, contradicting clinical and experimental observation for tumors.In contrast, by incorporating quiescent cells, we correctly predict that the core of the tumor is necrotic, that the outer layer of cells is fastest growing and most metabolically active, and that the tumor must reach an asymptotic size.
For the special case in which the transition probability from quiescent to proliferative cells is either zero or a constant value, representing a strict threshold or step function, our equations simplify.Over a given time interval, dt, we can integrate over the increased volume during that time, as determined by the thresholds calculated above, such that f p→q m p (t) = router,q(t+dt) router,q(t) f p→q (r)m p (r) dr, The probability density where B c (r) is the power intake of a cell, given a certain radius, and Q is a normalization constant.Because f is a probability density, its total probability over its support must equal 1.Therefore, Q router,q(t+dt) router,q(t) dr = 1 Thus, Q = 1 r outer,q (t + dt) − r outer,q (t) .
Using the fact that mass is the product of density and volume, we have router,q(t+dt) = 4 3 πρ T r outer,q (t + dt) 3 − r outer,q (t) 3 r outer,q (t + dt) − r outer,q (t) An exactly analogous formula holds for the necrosis rate of quiescent cells, Γ q , obtained by replacing r outer,q → r outer,n .Again, these growth dynamics differ fundamentally from those found in healthy tissue [97] where necrotic cores do not develop and cells rarely become quiescent due to lack of resources.Moreover, this formulation has the advantage of being more spatially explicit than the formulation given in Sec. ( 4) and [31].
6. Timing of drug therapies and effects of quiescent cells.New effects arise when drugs are used to treat tumors [49,50,51,52,53,54,88,73,2].In particular, we add two new terms to our coupled differential equations.One term represents direct death of proliferative cells due to drugs at a rate γ p .These cells apoptose and are cleared, so do not become part of the necrotic core, and thus exit from the equations altogether.This possibility is now allowed because drugs target and kill actively dividing cells.The other term allows for quiescent cells to return to a proliferative state at a rate f q→p .Without drugs this situation is extremely unlikely because once cells are below the threshold for available energy, it is unlikely they will ever gain access to more resources in the future.Drugs, however, will kill the outer layer of proliferative cells, allowing the quiescent cells to once again gain direct and easy access to blood supply and resources, thus returning to their proliferative state.
d m p (t) The timing and concentration of drug delivery enters through the parameter γ p , and we solve these equations numerically to explore three common examples.
Scenario 1: One pulse of the drug results in a concentration that decreases exponentially in time as d = d 0 e −νt .Eventually, the drug concentration will be low, and the tumor will be able to grow as it does in the radial dependence model.
Scenario 2: Multiple injections of the drug result in a saw function where decay is assumed to be exponential (as in Scenario 1) after the pulse of each injection.
Scenario 3: Slow-release drugs (idealized) that can be modeled either by a constant function or the saw function with a slower rate of decay.
In addition, we begin to incorporate the effects of pressure into our model [5,84].High pressures can squeeze shut blood vessels, preventing both the flow of oxygen and resources to fuel tumor growth as well as the delivery of drugs to kill dividing tumor cells.Therefore, there are two reservoirs that can allow tumors to be resilient and re-grow after treatment.First, quiescent cells can switch back to a proliferative state once drugs have eliminated proliferative cells that compete for resources and the drugs diminish to a non-effective concentration.Second, blood vessels closing due to high pressure, or simply patches of the tumor without blood vessels, can allow tumor cells to avoid the drugs altogether.
The full effects of pressure within the tumor are still being studied, but a useful function for modeling relative pressure was first derived and tested by Baxter and Jain [5].We modify their functional form to be in terms of absolute pressure is the surface area per unit volume for transport in the tumor, p o is outside pressure, and p n is a normalization factor for the pressure.According to this equation, the greater the pressure, the higher the likelihood that the cell will become quiescent or necrotic.Power thresholds for transition to quiescence and necrosis are given by B c,q = 0.001 J/day and B c,n = 0.0007 J/day respectively.transition probabilities are the same as that mentioned in Section 5, with the one difference being that the lower bound of integration is the radius at which Bc, q/Bc, n = 1/p(r).Parameter values are given in the figure caption.For Scenario 1 (single drug pulse with exponential decay), following treatment and as the drug concentration falls, the quiescent cells that survived can transition back to a proliferative state if they have access to sufficient resources, and thus, the tumor regrows.For Scenario 2 (multiple injections), however, the mass of quiescent and necrotic cells is eventually depleted because quiescent cells that transition to become proliferative are quickly killed by the next injection of drugs.Moreover, the necrotic mass decreases because the necrotic cells are slowly cleared by the body, Figure 11.Plots for resultant trajectories of quiescent cell numbers in response to the corresponding (by row from top to bottom) sequence of drug treatments in Fig. 9.These plots are model predictions based on numerical solutions obtained by evolving the equations through time.Parameter values are the same as are given in Fig. 10. and there are fewer and fewer quiescent cells available to transition to a necrotic state.For Scenario 3 (slow-release drugs), these effects result in a steady, monotonic decrease for all parts and cell states of the tumor.
Quiescent cells rarely transition back to a proliferative state in healthy tissue, and when it does occur, it often results in abnormal growths such as tumors.Consequently, this re-growth of the tumor after treatment with drugs is also a novel feature not typically observed or modeled for healthy tissue.Indeed, the growth equations for tumors represent novel and mathematically interesting features of growth beyond that usually expected for healthy tissue.That is, tumor growth represents an ideal case for applying more sophisticated mathematical tools to growth dynamics in biology.
7. Concluding remarks.In this paper we have developed new theory that explicitly connects fractal properties of tumor vasculature and growth with cell and host properties.In so doing, we can predict how tumor growth changes with tumor and host size, as well as in response to different timing sequences for treatment with drugs.In addition, we have contrasted the fractal characteristics of tumor vasculature and growth with healthy tissue and discovered intriguing differences.First, for tumor vasculature, the scaling exponents for vessel radii and lengths across branching levels are smaller, indicating a slower decrease in vessel dimensions, less hierarchy overall, and perhaps larger terminal vessels in tumors than in their healthy hosts.Second, and consistent with the first finding, despite the fact that tumor growth trajectories fall on the same universal curve as ontogenetic growth, tumors spend much more of their time in the exponential phase of growth and much less time in the asymptotic, logistic phase.Consequently, tumors may not make the transition to an asymptotic state that allows for the maintenance of existing cells.Indeed, tumors may grow as closely to exponentially as possible, maximizing the growth needs and fitness of each individual cell, until the host can no longer maintain this energetic burden and structural disturbance, eventually resulting in death.That is, by evading constraints, such as continued cell division and apoptosis, that are necessary for multicellularity, the tumor cells divide as quickly as possible with the resources available until the system breaks down.A critical future challenge is to more deeply understand how this exponential growth in tumors is maintained, and perhaps how to induce a transition within tumors to a logistic, asymptotic growth stage similar to that for ontogenetic growth.
To understand how drugs can be used to treat tumors, we need to model how the drugs affect all of the cell states-proliferative, quiescent, and necrotic-within the tumor.We find that quiescent cells, which arise only when resources are insufficient (often due to poor vascularization) for division, play a particularly important role for tumor re-growth after treatment.This finding is expected because drugs specifically target the proliferative cells and because necrotic cells are of course already dead.The advantage of our approach is that we can use a mathematical model to derive explicit quantitative predictions that should help clinicians anticipate what combinations of cell states and which drug sequences will likely lead to the fewest number of quiescent cells and/or smallest chances for tumor re-growth.The eventual hope is to develop models that can be tailored to specific and known vasculature and/or tissue constraints.
Our findings show that tumors grow exponentially for as long as possible but eventually reach a maximum size at which available local resources in the host cannot fuel further growth.The transition to this steady state size appears to occur much more rapidly in tumors than the more gradual asymptotic transition observed for ontogenetic growth (Fig. 5).A complicating factor to these questions is the novel feature of metastasis for tumors.Metastasis is the process by which cells from a primary tumor enter the bloodstream and are transported to other parts of the body where they attempt to establish a new tumor [59,69].As for the primary tumor, these metastatic tumors need to develop their own blood supply and vascular system through angiogenesis to fuel the growth of these new, additional tumors.Metastasis is the primary distinction between benign and malignant tumors.It is also a process that enables the tumor to effectively circumvent its maximum possible size by establishing tumors in new locations in the host when local resources are depleted in its current location.Indeed, it is very likely that there are selective pressures on genotypes within the tumor to develop the ability to metastasize, and thus, metastasis is another reason to consider the use of models from evolution (selection on multiple genotypes) and ecology (competition, composition, and diversity of genotypes) to understand tumor growth dynamics for a range of genotypes, represented by different cell types and more than the three cell types presented here [19,64,40,66].These questions are beyond the scope of this paper but will likely be an essential part of future progress on understanding tumor growth and treatment.The importance of metastasis is also important in considering drug treatments because of the need to treat tumors during early stages of growth and before metastasis occurs.
In summary, mathematical models that tie together multiple aspects of tumor growth and treatment have a strong potential to yield quantitative insights that are of clinical significance.In the past century, researchers have witnessed a great deal of empirical work and progress on cancer.Pioneering studies [93,92,1,56,48,30] that use mathematical modeling to study tumors have laid the groundwork for current efforts at modeling the growth dynamics of tumors.At present, there is a great opportunity to use recent advances in mathematical modeling, as contained within the pages of this journal, to devise methods to control and counteract tumor growth.Indeed, some clinicians believe that controlling tumor growth will be a more attainable and thus successful strategy than current treatments for tumors.Here, we focused on how tumor vasculature, host properties, and timing of drug sequences affect the growth and treatment of tumors.Moreover, we explored how vasculature and growth in tumors compares with that in healthy tissues.Our results reveal compelling similarities and differences that suggest avenues for future research.
1078 V. M. SAVAGE, A. B. HERMAN, G. B. WEST AND K. LEU

Figure 2 .
Figure 2. Plot of the logarithm of the length of vessels, l k , versus the branching order, k, as determined by the Strahler method.Data are from[37,99,100,45,35,56,55] from 7 healthy sections (e.g., left anterior descending artery, pulmonary) of vascular systems in pig, cat, and human.For healthy tissue, the correlations are excellent and correspond to scaling exponents that on average are close to the value of b = 1/2 (see Fig.6), somewhat higher than predicted.See Fig.6for estimates of scaling exponents based on best fit regression lines as well as p-values, 95% confidence intervals, and r 2 values.

Figure 4 .
Figure 4. Plot of the logarithm of the radii of vessels, r k , versus the branching order, k, as determined by the Strahler method.Data are from[37,99,100,45,35,56,55] from 7 healthy sections (e.g., left anterior descending artery, pulmonary) vascular systems in pig, cat, and human.For healthy tissue, the correlations are excellent and correspond to scaling exponents that on average are close to the predicted value of a = 1/2 for larger vessels (see Fig.6).Also see Fig.6for estimates of scaling exponents based on best fit regression lines as well as p-values, 95% confidence intervals, and r 2 values.

Figure 5 .
Figure5.Plot of the logarithm of the radii of vessels, r k , versus the branching order, k, as determined by the Strahler method.Data are from tumors (mammary and colorectal carcinomas (A91 and A97)) in rat and human.For tumors, the correlations are weaker and the scaling exponents are significantly shallower than for healthy tissue shown in Fig.4.For both healthy tissue and tumors, a self-similar fractal network is a reasonable approximation for the vascular systems, but for tumors there still remains significant variation to explain.See Fig.6for estimates of scaling exponents based on best fit regression lines as well as p-values, 95% confidence intervals, and r 2 values.

Figure 6 .
Figure 6.Table of calculated scaling parameters and correlation coefficients for vascular systems based on data in Figs. 1 − 5. Average values for healthy tissue and tumors are shown in bold.Tumor vasculature exhibits higher variation (with p-values for length scaling exponents all being > 0.0001), possibly corresponding to noise in gene expression or chromosome segregation, and also shallower scaling exponents, representing slower scaling and possible representing a breakdown in vessel hierarchy.

Figure 7 .
Figure 7. Plots of dimensionless ratio versus dimensionless time as defined by Eqs.(29)-(30).Data in grey are for ontogenetic growth from 13 species of animals (see[98] for original data sources), ranging from guppy to cod to guinea pig, and data in color with square symbols are for tumor growth trajectories for C3H mammary carcinoma (dark blue), EMT6 mammary carcinoma (dark red), KHJJ mammary carcinoma (light green), NCTC (dark green), Flank (yellow), Primary fibroadenoma (red), Primary Osteosarcoma (light blue), and Walker Carcinoma (purple) tumor types (see[31,92] for original data references).This plot shows the extent to which all the growth data can be described by a single universal curve[98,24], and thus, the extent to which tumor growth is really analogous to ontogenetic growth form birth to adult.Fitted parameter values are taken from the original references.
p is the hydraulic conductivity coefficient ( cm mmHG * s ) of the microvascular wall, L i is the hydraulic conductivity ( cm 2 mmHG * s ) of the interstitium, S V

Figure 9 .Figures 9 −
Figure 9. Timing sequences for drug delivery for the three scenarios described in the text.Parameters for the drugs are d = 100 units per volume, threshold for cell death is 10 units per volume, and ∆t = 7x10 −4 days is the time step.

Figure 10 .
Figure 10.Plots for resultant trajectories of proliferative cell numbers in response to the corresponding (by row from top to bottom) sequence of drug treatments in Fig. 9.These plots are model predictions based on numerical solutions obtained by evolving the equations through time.Parameter values are A = B 0 m c /E c = 0.305 (1/day), β = 1 2a+b = 1, F = B m /E c = 0.14 (1/day), λ = 0, α = 1.4,and ν = 0.01 (1/day).Parameters for pressure are, κ = 36.823(taken from [5]), p n = 10 mm Hg, and p o = 3mm Hg.Power thresholds for transition to quiescence and necrosis are given by B c,q = 0.001 J/day and B c,n = 0.0007 J/day respectively.

Figure 12 .
Figure 12.Plots for resultant trajectories of necrotic cell numbers in response to the corresponding (by row from top to bottom) sequence of drug treatments in Fig. 9.These plots are model predictions based on numerical solutions obtained by evolving the equations through time.Parameter values are the same as are given in Fig. 10.