Tumor growth under angiogenic signaling – data fitting (part 1)

In todays post I will show how to fit the model to the experimental data using MATLAB with Optimization Toolbox. It will be one of the many posts in which I will try to illustrate problems with parameter estimation. Let’s start with introducing the model.

Hahnfeldt and colleagues proposed a mathematical model of the concept that tumor growth and host blood vessel support is bidirectionally modulated (Hahnfeldt et al., Cancer Research, 1999). Tumor volume (V) and effective vascular support (K) that defines the tumor carrying capacity are time-dependent variables described by a set of coupled ordinary differential equations (ODEs). Tumor growth is assumed to be governed by the Gompertz law:

With constant effective vascular support K=K_max, initial rapid tumor growth is followed by a slowdown as the tumor volume approaches carrying capacity K_max. To account for reciprocal interaction of the tumor with the host vasculature, carrying capacity through vascular support can be described by a variable K modulated by the tumor V:

where  represents spontaneous loss of functional vasculature,  represents vessels growth stimulation due to factors secreted by the tumor proportionally to its size, and describes endogenous inhibition of previously generated vasculature due to factors secreted by the tumor proportionally to the tumor surface-to-volume ratio.

We will try to see how well the model works on the actual experimental data (not the one from the paper). Quick Google image search using the term “mouse bevacizumab” (bavacizumab is antiangiogenic drug) returned a figure from the paper by Tsukihara an others (Tsukihara et al., Anticancer Research, 2015) in which we have a plot of relative tumor volume in time in the presence or absence of bevacizumab. We will use that figure to grab the experimental data. There is a very good tool available on MATLAB Central called “grabit” that allows to easily collect the data from the plot (screenshot below).

After calibrating the the axis, we can grab the average values and the standard deviation (for each point first grab the average then the top of the standard deviation bar, skip the point at day 0) first for the control experiment and than for Bevacizumab alone. After that it is enough to save the data into the file e.g. measurements.mat. First we need to write function that will read saved file and restructure it to the format that we want to work on.

function data = processMeasurements( file )
    %load the file with the measurements
    tmp = load(file); 
    %x values are in the first column and we need to take 
    %every second point (we have taken average and standard deviation at the same point).
    %We also round the values to the full days
    data.t = round(tmp.measurements(1:2:end/2,1))'; 
    %now we grab the averages. In first row we will have data for control
    %and in the second row for bevacizumab
    data.average = reshape(tmp.measurements(1:2:end,2),[],2)';
    %now we grab the standard deviations. In first row we will have data for control
    %and in the second row for bevacizumab
    data.std = reshape(tmp.measurements(2:2:end,2)-tmp.measurements(1:2:end,2),[],2)';
end

To make sure that everything is structured properly we can write the function that will plot grabbed data,

function plotData(data)
    figure(1)
    clf
    hold on
    errorbar(data.t,data.average(1,:), data.std(1,:),'Color','r','LineStyle','none','Marker','o');
    errorbar(data.t,data.average(2,:), data.std(2,:),'Color','b','LineStyle','none','Marker','s');
    hold off
    xlabel('day');
    ylabel('RTV (mean and std)');
    legend({'Control','Bevacizumab'},'Location','NorthWest')
end

which should result in the following plot

Now we can start working on the code that will fit the Hahnfeldt model to the data. In the part 1 of the post we will focus only on fitting to the control data – fitting both curves will be covered in part 2. The main script that we will run is the following

%read and restructure experimental data
data = processMeasurements('measurements.mat');

%set initial values for parameters and initial conditions
params = initializeParams();

%set which params should be estimated
%names need to be consistent with the params structure!
%in most of the studies based on Hahnfeldt model lambda2 is assumed to be 0
%we will do the same and won't fit its value
paramsToFit = {'lambda1','b','d'};

%perform data fitting
[params, err] = fitParameters(params,data,paramsToFit);

%solve the model for new parameters
sol = solveModel(data.t(end),params);

%plot data together with solution
plotDataAndSol(data, sol);

We have already implemented processMeasurements function, so we can write down the initializeParams function.

function params = initializeParams()

    %model parameters (we take the values from Hahnfeldt et al. paper)
    params.lambda1 = 0.192;
    params.lambda2 = 0;
    params.b = 5.85;
    params.d = 0.00873;
    
    %initialConditions
    %in the experimental paper implanted tumors had about 8 cubic milimiters
    params.V0 = 8;
    %initial tumor is not visualized. It is typically assumed that avascular tumor can grow to about 1 mm in radius
    params.K0 = 4/3*pi;
    
end

Now we need to code the fitParameters function. We will use the lsqnonlin function that is provided in Optimization Toolbox, as it handles by default constrains on the parameters values (we assume that all parameters are non-negative).

function [params, err] = fitParameters(params,data,paramsToFit)

    %preparing initial parameters values
    x0 = zeros(size(paramsToFit));
    for i = 1:length(paramsToFit)
        x0(i) = params.(paramsToFit{i});
    end
    
    opt = optimset('Display','iter');
    %Fmin is the function that returns the fit error for a given set of parameters
    %it will need data, params and paramsToFit variables so we pass them
    %as additional arguments to the lsqnonlin function 
    %(they will be passed forward to Fmin)
    flag = 0;
    while flag<1 %repeat until solver reached tolerance or converged
        opt = optimset('Display','iter');
        [x, fval, ~, flag] = lsqnonlin(@Fmin,x0,zeros(size(x0,[],opt,data,params,paramsToFit);
        x0 = x;
    end
    
    %rewritng x (vector of fitted values) back into params structure
    for i = 1:length(paramsToFit)
        params.(paramsToFit{i}) = x(i);
    end
    
end

The Fmin function will calculate the sum of squared differences between the model solution and the experimental data.

function err = Fmin( x, data, params,paramsToFit )

    %rewritng x into params structure
    for i = 1:length(paramsToFit)
        params.(paramsToFit{i}) = x(i);
    end
  
        %we solve the model
        sol = solveModel([0 data.t],params);
        %we calculate relative tumor volume at points of measurements
        RTV = sol.y(1,2:end)./sol.y(1,1);
        %we calculate the error
        err = RTV - data.average(1,:);

end

Now we need to write the solveModel function that returns model solution for given set of parameters and initial conditions.

function sol = solveModel( T, params )
    %we assume that initial conditions are in the params structure

    %we use MATLAB built-in ODE solver to calculate the solution
    sol = ode45(@ODE,[0 T(end)],[params.V0 params.K0]);

    %if user provided time mesh instead of only final point
    %we evaluate the solution on that mesh
    if length(T)>1
        sol.y = deval(sol,T);
        sol.x = T;
    end

    %model definition
    function dy = ODE(~,y)
       dy = zeros(2,1);
       dy(1) = -params.lambda1*y(1)*log(y(1)/y(2));
       dy(2) = -params.lambda2*y(2)+params.b*y(1)-params.d*y(2)*y(1)^(2/3);
    end

end

Finally, the function that will plot the solution and the data on the same plot.

function plotDataAndSol(data,sol)

    figure(2)
    clf
    hold on
    errorbar(data.t,data.average(1,:), data.std(1,:),'Color','r','LineStyle','none','Marker','o');
    plot(sol.x,sol.y(1,:)./sol.y(1,1),'k');
    hold off
    xlabel('day');
    ylabel('RTV (mean and std)');
    legend({'Data','Model'})

end

Now we are ready to run the main script. This results in the following plot, which shows nice correspondence of the model solution to the experimental data,

with estimated parameters values . It seems that there is nothing more to do – we estimated the parameters and nicely fitted model to the experimental data. However…

Let’s assume that we didn’t grab the data from the plot – we don’t consider that source of the noise. We all know, however, that there is uncertainty in experimental measurements, which is hard to estimate. Let us perform the same fitting procedure, but for the data that is perturbed by up to 1% (each point, uniformly) – quite a small measurement error!

Plotting the estimated values of parameters (see below) shows that even such a small error introduces huge variability with estimated parameters values. Parameter b has the values in the range from 384 to 1720 and d in the range from 3 to 11. How certain can we be with our parameters estimations then? That is quite a simple model and “clean” experimental data…

Advertisement