Capturing the Driving Role of Tumor-host Crosstalk in a Dynamical Model of Tumor Growth

下载 PDF 引用 收藏 提问与回复 分享您的反馈 Cited by



Cancer Research
Mar 2015


In 1999, Hahnfeldt et al. proposed a mathematical model for tumor growth as dictated by reciprocal communications between tumor and its associated vasculature, introducing the idea that a tumor is supported by a dynamic, rather than a static, carrying capacity. In this original paper, the carrying capacity was equated with the variable tumor vascular support resulting from the net effect of tumor-derived angiogenesis stimulators and inhibitors. This dynamic carrying capacity model was further abstracted and developed in our recent publication to depict the more general situation where there is an interaction between the tumor and its supportive host tissue; in that case, as a function of host aging (Benzekry et al., 2014). This allowed us to predict a range of host changes that may be occurring with age that impact tumor dynamics. More generally, the basic formalism described here can be (and has been), extended to the therapeutic context using additional optimization criteria (Hahnfeldt et al., 1999). The model depends on three parameters: One for the tumor cell proliferation kinetics, one for the stimulation of the stromal support, and one for its inhibition, as well as two initial conditions. We describe here the numerical method to estimate these parameters from longitudinal tumor volume measurements.

Keywords: Mathematical oncology (数学科), Tumor growth (肿瘤的生长), Parameters estimation (参数估计)


  1. Matlab (The Mathworks Inc) version 12 or later with optimization toolbox (can also be adapted for Scilab: http://www.scilab.org/index.php/en)


This protocol describes how to use a dedicated mathematical model for analysis of tumor growth kinetics. Specifically, it deals with a method for determination of the coefficients of the model from longitudinal measurements of tumor growth. In the data example provided, the measures were obtained by using calipers to determine, every other day, largest (L) and smallest (w) diameters of subcutaneously implanted tumors (Benzekry et al., 2014) for additional details. The formula V = w2Lπ/6 was then used to compute the tumor volume. Researchers can use this protocol in order to compare two groups (or more) of tumor kinetics data and identify which of the coefficient(s) of the mathematical model significantly differ(s) (or not) among the groups. Longitudinal measurements of growth data are required.
A two-dimensional ordinary differential equation model for tumor growth as a function of the tumor volume, V, and its carrying capacity, K (formally, the maximum tumor volume that level of host environmental support could sustain), was used (t denotes time):

with V0 = the volume of the tumor at t = t0, K0= the carrying capacity at t = t0, a = proliferation of the tumor cells, b = carrying capacity stimulation, and d = carrying capacity inhibition. Numerical solutions of the model are fitted to experimental measurements of tumor volume at several time points, and the set of parameters (a, b, d) generating the best fit of the model are thereby determined. In other words, the values of (a, b, d) are determined as the ones minimizing the distance between the model simulation and the data (thus maximizing the likelihood of the data under the hypothesis that it has been generated by the model). The quantity V0 is fixed to the first observed data point and K0 is fixed to 2V0. For instance, for the first animal of group 1 of the data example, the resulting values of the parameters were: V0 = 836 mm3, K0 = 1672 mm3, a = 0.224 day-1, b = 0.710 day-1 and d = 0.0018 mm-2 day-1.
This protocol is composed of several Matlab scripts and functions that are freely downloadable at the following address: https://github.com/benzekry/fit_tumor_growth. In the following, we detail each step that should be sequentially launched (by typing the script name in a Matlab command window) in order to perform the full task. The associated script is indicated in bold font. We perform these on an illustrative data set example composed of two tumor growth curve groups, each n = 20 mice.

  1. Import the data from Excel (xlsx) file into Matlab. importData.m.
  2. Manually explore the parameter space to determine an initial parameter guess that approximately fits the data. play_model.m.
  3. (Optional) Determine bounds on the parameters. If the parameter space is not restricted for this model, one may face identifiability issues (several parameter sets yielding almost equal solutions). In our example, we chose bounds so that each parameter spans two orders of magnitude around the initial guess. This is specified in fitGlobal.m (within the function lsqcurvefit).
  4. Determine the error variance model to be employed (uncertainty on the data). Options include proportional (common) or constant (unadvised for tumor growth measurements). Note that the former is equivalent to a fit (i.e., sum of squared residuals minimization) of the log-transformed data by the log-transformed model. For caliper-derived measurements of tumor growth, a slightly more elaborate error model might be employed (Ledzewicz and Schattler, 2008). It involves a power α < 1 of the data as well as a minimal threshold and leads to weighted least-squares. In our example we used a proportional error model (specified in lsqcurvefit within fitGlobal.m).
  5. Fit the model to the data using a built-in minimization function of Matlab. Options incude: lsqcurvefit, fmincon, and fminsearch. Personal experience suggests a better ability of fminsearch to escape local minima, although no bound constraint on the values of the parameters can be prescribed with this function. This can be circumvented by penalizing the model for parameter values outside the bounds. Alternatively, stochastic algorithms (such as Markov Chain Monte Carlo) or multi-start optimization might also be considered. We recommend fitting each animal growth curve separately, in particular when inter-animal variability is high, as opposed to fitting the average or median growth curve. The script that launches the fit and exports individual plots is launch_fit.m. Individual plots are exported into folders that have user-defined names (‘Group1’ and ‘Group2’ in our example).
  6. Check the goodness-of-fit. This includes:
    1. Visual inspection of the simulated individual growth curves against the data (Figure 1).
    2. Distribution of the (weighted) residuals (all individuals pooled together). It should be gaussian, a hypothesis that can be assessed using the Kolmogorov-Smirnov test (as can be performed using the Matlab function kstest.m) (Massey, 1951). Two plots of the residuals are outputs when launching launch_fit.m.
    3. Computation of the coefficient of determination (R2), defined by , where yi is the data point at time ti, f(ti) is the solution of the model at this time and is the average of the yi’s. The value of R2 is printed at the screen by launch_fit.m. This number is a usual metric of goodness-of-fit that quantifies how better is the model at explaining the data in comparison to the average only (which would give R2 = 0). The closer R2 is to one, the better the fit.

    Figure 1. Individual fits of the tumor growth data by the model for each mouse from one of the two groups. The blue curve is the volume (V) curve while the red curve depicts the dynamics of the carrying capacity (K). Fits were performed by optimization of the coefficients a, b and d. Initial volume V0 was set to the first volume measured. Initial carrying capacity was set to the double of this quantity. x-axis = time in days, y-axis = tumor size (in mm3).

  7. Investigate for statistical differences in the parameters of the model (here a, b, d and K0) and across the groups under study. File manips_postFitAnalysis.m outputs two things:
    1. Boxplots of the parameter values, including statistical significance (t-test) (Figure 2). Statistical power depends on the number of mice employed. With n = 20 animals, the analysis performed here was able to detect a significant difference of 0.205 day-1 in parameter a between the two groups (0.295 ± 0.0758 in group 1 versus 0.09 ± 0.0132 in group 2, mean ± standard error, see Figure 2), i.e., 69.4%, at the significance level of α = 0.05. To detect a similar difference of 0.2 day-1 with a probability of success (power) of 95% (i.e., probability of false positive of 5%), only 10 animals are sufficient (and necessary).

    Figure 2. Boxplots of the distribution of three parameters of the model (a, b and d) for the two groups. Initial condition V0 and carrying capacity K0 were fixed as in Figure 1. x-axis = Group 1 and Group 2. y-axis = parameter values for a, b and d. * = α < 0.05 (Student’s t-test).

    1. Simulated growth curves: Either all the animals or averaged other the groups (Figure 3).

    Figure 3. Simulated growth curves. Left panel: All individual simulations corresponding to each animal-specific fitted parameter set. Right panel: Average growth curves (mean ± standard error). x-axis = time in days, y-axis = tumor volume in mm3. Blue curves = Group 1. Red curves = Group 2.


This project was supported, in part, by the National Aeronautics and Space Administration under NSCOR grants NNJ06HA28G and NNX11AK26G and by Award Number U54CA149233 from the National Cancer Institute, both to L. Hlatky. This study also received support within the frame of the LABEX TRAIL, ANR-10-LABX-0057 with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the “Investments for the future” Programme IdEx (ANR 10-IDEX-03-02).


  1. Beheshti, A., Benzekry, S., McDonald, J. T., Ma, L., Peluso, M., Hahnfeldt, P. and Hlatky, L. (2015). Host age is a systemic regulator of gene expression impacting cancer progression. Cancer Res 75(6): 1134-1143.
  2. Benzekry, S., Lamont, C., Beheshti, A., Tracz, A., Ebos, J. M., Hlatky, L. and Hahnfeldt, P. (2014). Classical mathematical models for description and prediction of experimental tumor growth. PLoS Comput Biol 10(8): e1003800.
  3. Hahnfeldt, P., Panigrahy, D., Folkman, J. and Hlatky, L. (1999). Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy. Cancer Res 59(19): 4770-4775.
  4. Ledzewicz, U. and Schattler, H. (2008). Optimal and suboptimal protocols for a class of mathematical models of tumor anti-angiogenesis. J Theor Biol 252(2): 295-312.
  5. Massey, F. J., Jr. (1951). The Kolmogorov-Smirnov test for goodness of fit. J Am Statist Assoc 46(253): 68-78.


1999年,Hahnfeldt等人提出了肿瘤生长的数学模型,其由肿瘤及其相关脉管系统之间的相互通信决定,引入了这样的想法:肿瘤由动态而非静态的支持,承载能力。在该原始文件中,携带能力等同于源自肿瘤衍生的血管生成刺激剂和抑制剂的净效应的可变肿瘤血管支持。这种动态承载能力模型在我们最近的出版物中被进一步抽象和开发,以描述在肿瘤与其支持性宿主组织之间存在相互作用的更一般的情况;在这种情况下,作为宿主老化的函数(Benzekry et al。,2014)。这使我们能够预测可能随着年龄发生的影响肿瘤动力学的宿主变化的范围。更一般地,这里描述的基本形式可以(并且已经)使用附加的优化标准扩展到治疗上下文(Hahnfeldt等人,1999)。该模型取决于三个参数:一个用于肿瘤细胞增殖动力学,一个用于刺激基质支持,一个用于其抑制,以及两个初始条件。我们在这里描述数值方法从纵向肿瘤体积测量估计这些参数。

关键字:数学科, 肿瘤的生长, 参数估计


  1. Matlab(The Mathworks Inc)版本12或更高版本与优化工具箱(也可以适应Scilab: http: //www.scilab.org/index.php/en


该协议描述如何使用专用的数学模型来分析肿瘤生长动力学。具体地,其涉及用于从肿瘤生长的纵向测量确定模型的系数的方法。在所提供的数据示例中,通过使用测径器每隔一天确定皮下植入的肿瘤的最大(em)和最小(em)w直径(Benzekry et al。,2014)了解更多详情。然后使用公式V = w 2 Lπ/6 来计算肿瘤体积。研究人员可以使用该方案以便比较两组(或更多)肿瘤动力学数据并且识别数学模型的系数中的哪一个在组之间显着不同(或不同)。需要纵向测量生长数据 作为肿瘤体积的函数的肿瘤生长的二维常微分方程模型,和其携带能力,K em(形式上,肿瘤体积的最大肿瘤体积的主机环境支持可以维持),被使用( t 表示时间):

与 t 的肿瘤体积 > 0 = t a =肿瘤细胞增殖,能力刺激,以及 d =携带能力抑制。模型的数值解决方案适合于在几个时间点的肿瘤体积的实验测量,并且参数集合(a,b,d,d,d, )产生模型的最佳拟合。换句话说,( a , b , d 的值被确定为最小化模型模拟和数据之间的距离(从而使数据在由模型生成的假设下的可能性最大化)。数量 V 0 固定为第一个观测数据点和 K /em> 固定为 2 。例如,对于数据示例的组1的第一动物,参数的结果值为: 3 , K 0 = 1672 mm 3 = 0.224天 -1 , b = 0.710天 -1 和 d = 0.0018mm 2 day -1
此协议由几个Matlab脚本和函数组成,可在以下地址免费下载: https://github.com/benzekry/fit_tumor_growth。在下面,我们详细说明每个应该顺序启动的步骤(通过在Matlab命令窗口中键入脚本名称),以执行完整的任务。相关的脚本用粗体字表示。我们对由两个肿瘤生长曲线组组成的说明性数据集示例执行这些,每个n = 20只小鼠。

  1. 将数据从Excel(xlsx)文件导入 Matlab。 importData.m 。
  2. 手动探索参数空间以确定近似拟合数据的初始参数猜测。 play_model.m 。
  3. (可选)确定参数的边界。如果参数空间不限于此模型,可能面临识别问题(几个参数集产生几乎相等的解)。在我们的例子中,我们选择了边界,使得每个参数围绕初始猜测跨越两个数量级。这在fitGlobal.m中指定(在函数 lsqcurvefit 内)。
  4. 确定要采用的误差方差模型(数据的不确定性)。选项包括比例(共同)或常数(未设定肿瘤生长测量)。注意,前者等价于对数变换模型的对数变换数据的拟合(即,即平方残差最小化的和)。对于肿瘤生长的卡尺导出的测量,可以采用稍微更复杂的误差模型(Ledzewicz和Schattler,2008)。它涉及功率α 1的数据以及最小阈值,并导致加权最小二乘法。在我们的示例中,我们使用了一个比例误差模型(在 fitGlobal.m 中的 lsqcurvefit 中指定)。
  5. 使用Matlab的内置最小化函数将模型拟合到数据。选项incude:lsqcurvefit,fmincon和fminsearch。个人经验表明fminsearch逃避局部最小值的更好的能力,虽然对于该函数不能规定参数值的约束约束。这可以通过对边界外的参数值的模型进行惩罚来规避。或者,也可以考虑随机算法(例如马尔可夫链蒙特卡罗)或多启动优化。我们建议分别拟合每个动物生长曲线,特别是当动物间变异性高时,而不是拟合平均或中间生长曲线。启动拟合和导出单个地块的脚本是 launch_fit.m 。单个图形导出到具有用户定义名称(本示例中为"Group1"和"Group2")的文件夹中。
  6. 检查拟合优度。这包括:
    1. 目视检查模拟个体生长曲线与数据(图1)。
    2. (加权)残差的分布(所有个体合并 一起)。它应该是高斯的,一个可以评估的假设 使用Kolmogorov-Smirnov检验(如可以使用Matlab进行的) 函数 kstest.m )(Massey,1951)。残差的两个图是 启动launch_fit.m时输出。
    3. 计算 由 2 " height ="50" width ="127"/>,其中 y i i sub> )是此模型的解决方案,是 y i 的平均值。 launch_fit.m 会在屏幕上打印 R 2 的值。这个数字是适合度的通常度量 量化了在比较解释数据时模型有多好 只有平均值(这将给予 = 0 )。更接近的 R2 是一个, 配合越好。

    图1.通过模型对来自两组之一的每只小鼠的肿瘤生长数据的个体拟合。蓝色曲线是体积( V )曲线,曲线描述了承载能力(K em)的动力学。通过优化系数a , b 和 d 进行拟合。将初始体积V 0设置为所测量的第一体积。初始承载力设定为该数量的两倍。 x轴=以天计的时间,y轴=肿瘤大小(以mm为单位)。

  7. 调查模型参数的统计差异(此处 a , b , d 和 K 0 )。文件manips_postFitAnalysis.m输出两个东西:
    1. 参数值的箱线图,包括统计显着性 (t检验)(图2)。统计力量取决于小鼠的数量 采用。对于n = 20只动物,在此进行的分析能够 检测参数a之间的0.205天 -1 的显着差异 两组(组1中为0.295±0.0758,组中为0.09±0.0132) 2,平均值±标准误差,参见图2),即,69.4%的显着性 ?α的水平<0.05>。为了检测0.2天 -1 与a的相似差异 成功概率(功率)为95%(即,假的概率) 阳性5%),只有10只动物足够(和必要)。

    图2.两组的模型(a,b和d)的三个参数的分布的箱线图。初始条件 V /em> 和承载能力 K 0 如图1所示。x轴=组1和组2。 y轴=
    , b 和 d 的参数值。 * =α< 0.05(学生t检验)。

    1. 模拟生长曲线:所有动物或平均其他组(图3)。

    图3.模拟生长曲线。左图:所有单个模拟对应于每个动物特异性拟合的参数集。右图:平均生长曲线(平均值±标准误差)。 x轴=以天计的时间,y轴=以mm 3为单位的肿瘤体积。蓝色曲线=组1.红色曲线=组2.


该项目得到国家航空航天局NSCOR授予NNJ06HA28G和NNX11AK26G的部分支持,以及来自国家癌症研究所的授予U54CA149233,均授予L. Hlatky。这项研究还得到了法国国家资助的LABEX TRAIL,ANR-10-LABX-0057框架内的支持,由法国国家研究机构(ANR)在"未来投资"计划框架下管理IdEx(ANR 10-IDEX-03-02)。


  1. Beheshti,A.,Benzekry,S.,McDonald,J.T.,Ma,L.,Peluso,M.,Hahnfeldt,P.and Hlatky, 宿主年龄是影响癌症进展的基因表达的系统调节因子 癌症Res 75(6):1134-1143。
  2. Benzekry,S.,Lamont,C.,Beheshti,A.,Tracz,A.,Ebos,J.M.,Hlatky,L。和Hahnfeldt,P.(2014)。 用于描述和预测实验性肿瘤生长的经典数学模型 PLoS计算Biol 10(8):e1003800。
  3. Hahnfeldt,P.,Panigrahy,D.,Folkman,J。和Hlatky,L。(1999)。 在血管生成信号传导下的肿瘤发展:肿瘤生长,治疗反应和血管后休眠的动力学理论。/a> Cancer Res 59(19):4770-4775
  4. Ledzewicz,U。和Schattler,H。(2008)。 针对肿瘤抗血管生成的一类数学模型的最优和次优方案。 em> J Theor Biol 252(2):295-312。
  5. Massey,F.J.,Jr.(1951)。 Kolmogorov-Smirnov测试的适合度。 J Am Statist Assoc 46(253):68-78。
  • English
  • 中文翻译
免责声明 × 为了向广大用户提供经翻译的内容,www.bio-protocol.org 采用人工翻译与计算机翻译结合的技术翻译了本文章。基于计算机的翻译质量再高,也不及 100% 的人工翻译的质量。为此,我们始终建议用户参考原始英文版本。 Bio-protocol., LLC对翻译版本的准确性不承担任何责任。
Copyright: © 2015 The Authors; exclusive licensee Bio-protocol LLC.
引用:Benzekry, S., Beheshti, A., Hahnfeldt, P. and Hlatky, L. (2015). Capturing the Driving Role of Tumor-host Crosstalk in a Dynamical Model of Tumor Growth. Bio-protocol 5(21): e1644. DOI: 10.21769/BioProtoc.1644.