sbiofit ability to estimate the parameter for 3 different cases at the same time

Hi,

I wrote the code below for estimation of KD from the binding assay for three different cell lines: i.e. bxpc3, colo205, T84 which have different number of antigen molecule on their surface. The code properly work for estimation of kD spearately for each cell line.

How can I estimate the kD as one value for different cell lines at the same time?

I attached my excel file (exp.xlsx) in which each cell line data has been stored in a spearted sheet.

Thanks

%% Setup Model
case_study = 'T84'; % 'bxpc3' 'colo205' 'T84'
model_input_file='exp.xlsx'
switch case_study   
    case 'bxpc3'
        exp_data=readtable(model_input_file,'Sheet','bxpc3');
        data.dose = exp_data.dose;
        data.ro = exp_data.RO_exp;
        input1=[0.1 0.4 1.2 11.1 33.3 100 300];
        RPC_Ag_Tumor = 219202;
    case 'colo205'
        exp_data=readtable(model_input_file,'Sheet','colo205');
        data.dose = exp_data.dose;
        data.ro = exp_data.RO_exp;
        input1=[0.1 0.4 1.2 3.7 11.1 33.3 100];
        RPC_Ag_Tumor = 662642;
    case 'T84'
        exp_data=readtable(model_input_file,'Sheet','T84');
        data.dose = exp_data.dose;
        data.ro = exp_data.RO_exp;   
        input1=[0.1 0.4 1.2 3.7 11.1 33.3 100 300];    
        RPC_Ag_Tumor = 169163;
end
data = structfun( @rmmissing , data , 'UniformOutput' , false);
setup; % (I build the model here)
%% Data Import
gData = groupedData(exp_data);
gData.Properties.VariableUnits = {'','second','nanomole','dimensionless'};
gData.Properties.GroupVariableName       = 'ID';
gData.Properties.IndependentVariableName = 'Time';
sbiotrellis(gData,'ID','Time',{'RO_exp'},'Marker','+','LineStyle','none');
% Extract data columns into separate variables for the weight expression
% evaluation.
headings = exp_data.Properties.VariableNames;
for i=1:length(headings)
    next = headings{i};
    eval([next ' = exp_data. ' next ';']);
end
%% Create the data dose.
d2                  = sbiodose('dose');
d2.TargetName       = 'Ab';
d2.LagParameterName = '';
d2                  = createDoses(gData, 'dose', '', d2);
% Build table of doses.
d2   = num2cell(d2);
dosesForFit = d2;
%% Fitting options
% Define response information.
responseMap = {'RO = RO_exp'};
% Define objects being estimated and their initial estimates.
estimatedInfoObj = estimatedInfo({'kD_Ab'});
% estimatedInfoObj = estimatedInfo({'k_max', 'EC50'});

estimatedInfoObj(1).Bounds       = [1E-1 1E2];
% estimatedInfoObj(2).Bounds       = [1E-6 1E-4];

% Define Algorithm options.
options                     = optimoptions('particleswarm');
options.Display= 'iter'; 
% Define fit problem.
f              = fitproblem('FitFunction', 'sbiofit');
f.Data         = gData;
f.Model        = model;
f.Estimated    = estimatedInfoObj;
f.ErrorModel = 'exponential';
f.ResponseMap  = responseMap;
f.Weights      = [];
f.Doses        = dosesForFit;
f.FunctionName = 'particleswarm';
f.Options      = options;
f.ProgressPlot = true;
f.UseParallel  = false;
f.Pooled       = true;
%% Run the model with optimum parameters
fitResults = f.fit;
plot(fitResults)
ci = sbiopredictionci(fitResults);
plot(ci)
ciParam = sbioparameterci(fitResults);
ci2table(ciParam)
plot(ciParam)



fitResults.ParameterEstimates
%update parameters
for i=1:length(estimatedInfoObj)
    xopt = sbioselect(model,'Name',fitResults.EstimatedParameterNames{i});
    xopt.value = fitResults.ParameterEstimates.Estimate(i);
end
params = {'Ab0'};
obs = {'RO'};
input=dose(~isnan(dose));
RO_exp = RO_exp(~isnan(RO_exp));
sfxn = createSimFunction(model, params, obs, d1,'UseParallel',false,'AutoAccelerate',false);
sfxn.accelerate;
doseTable = getTable(d1);
simData=sfxn(input,configsetObj.StopTime,doseTable);
%% Plot
fontsize = 18;
linesize = 2;
markersize = 10;
f = figure('Position', [14 10 900 600]);
colors = ["k", "#D95319", "#A2142F", "#77AC30", "#7E2F8E", "#00CED1", "#FF1493"];  % Define an array of colors
set(f, 'Visible', 'on');
set(gca,'fontsize', fontsize); 
hold on;
box on;
grid on;
set(gca, 'XScale', 'log');
%set(gca, 'YScale', 'log');
for i = 1:size(input,1)
   if i < size(input,1)
       line([input1(i) input1(i+1)], [simData(i).Data(end) simData(i+1).Data(end)], 'Color', colors(3),'LineWidth', 2);
   end
   % plot(input1(i), simData(i).Data(end), 'O', 'MarkerEdgeColor', colors(3), 'MarkerFaceColor', colors(3), 'markersize', markersize);

   plot(input1(i), RO_exp(i), 'O', 'MarkerEdgeColor', colors(4), 'MarkerFaceColor', colors(4), 'markersize', markersize);
end
xlabel('Dose (ug/mL)', 'fontweight', 'bold', 'Fontsize', fontsize, 'Interpreter', 'none');
ylabel('RO (%)', 'fontweight', 'bold', 'Fontsize', fontsize);
legend('Simulation','Experiment' , 'NumColumns', 1, 'Location', 'best', 'FontSize', fontsize)
saveas(gca, fullfile(save_locations{2}, 'RO'), 'png');
close(f);

Answer 1

[Jeremy Huard]

Hi @Mehdi,

Currently, when you run this code, you read one sheet and save the experimental data into the variable exp_data, which is then converted into a groupData object and used for fitting.

If you want to use all sheets, you will need to read them all and concatenate the resulting tables vertically to get one single data set.

If the dose targets across cell lines are the same, you can use the same ID for all experiments with the same dose amount.

And because you are already using the option Pooled=true, one single value for kD will be estimated across all groups.

Best regards,

Jérémy

[Mehdi]

Thanks for reply. Considering the structure of exp.xlsx in which some rows does not have any value and also different values of antigen experssion for each cell line RPC_Ag_Tumor and input1, can you suggest how I can implement your your solution?

Thanks

[Jeremy Huard]

I don't see where RPC_Ag_Tumor and input1 are used. RPC_Ag_Tumor is only defined but does not seem to be used for fitting. input1 is also only used for plotting, not for fitting.

That said, input1 seems to be the dose contained in the Excel sheet converted into some other units. Is this correct?

[Mehdi]

As you can see, once I call 'seup' I am building the model as below:

setup subroutine:

%% Model Name and Version
v = string(1); %version
model_name = 'binding'+v;
model_input_file='exp.xlsx';
%% Model Settings
configurations;
%% Parameters
parameters;
%% Compartments
compartments;
%% Species
species;
%% Rules
rules;
%% Reactions
reactions;
%% Dose Schedule
dose_schedule;
%% Verification
verify (model);
%% SBML and Excel Output
export;

in parametersubroutine in setup the RPC_Ag_Tumor is used as follows:

Avo = 6.022E14; % Avogadro number*1E-9
paramInfo = {
    'kon_Ab', 5E-4,  'liter/nanomole/second',  1, 'Association rate constant for binding of ADC or Ab with Ag';
    'kD_Ab', 1,  'nanomole/liter',  1, 'Equilibrium binding constant for Ag:Ab complex';
    'koff_Ab', 1E-5,  '1/second',  1, 'Initial Assignment';
    'kdeg', 1E-4,  '1/second',  1, 'Degradation rate constant for Ag';
    'Vcell', 1E-12,  'liter',  1, 'Volume of a cell';
    'Ncell', 1E7,  'dimensionless',  1, 'Number of cells in an in vitro experiment';
    'RPC_Ag_Tumor', RPC_Ag_Tumor/Avo,  'nanomole',  1, 'Receptor per cell expressed on a tumor cell';
    'Ag_total', 0,  'nanomole',  1, 'Total Initial Ag in well';
    'ksyn', 0,  'nanomole/second',  1, 'Initial Assignment';
    'Ab0', 0.005,'nanomole',1,'Ab initial concentration';
    'RO', 0,  'dimensionless',  0, 'Receptor Occupency';
    };
% Loop through the parameters and create the parameter objects
for i = 1:size(paramInfo, 1)
    % Create the parameter object using addparameter function
    parameterObj = addparameter(model, paramInfo{i, 1}, paramInfo{i, 2}, 'ValueUnits', paramInfo{i, 3},'ConstantValue', paramInfo{i, 4},'Notes', paramInfo{i, 5});
end	

but input1 is only for plotting, you are right.

[Jeremy Huard]

Hi @Mehdi,

I would use group-specific variants in that case where group-specific parameter values are defined in the dataset itself in a new column.

In your case, it could look like this:

% Import dataset into one single table and include group-specific values of RPC_Ag_Tumor

model_input_file='exp.xlsx';

snames = sheetnames(model_input_file);

Avo = 6.02214076e23;

RPC_Ag_Tumor = [219202, 662642, 169163] / Avo;

exp_data = table;

for sn=1:numel(snames)

exp_data_temp = readtable(model_input_file,Sheet=snames(sn));

exp_data_temp.RPC_Ag_Tumor(:) = RPC_Ag_Tumor(sn);

if sn >= 2

exp_data_temp.ID = exp_data_temp.ID + max(exp_data.ID);

end

exp_data = [exp_data; exp_data_temp]; %#ok<AGROW>

end

% model setup

setup;

% Data Import

gData = groupedData(exp_data);

gData.Properties.VariableUnits = {'','second','nanomole','dimensionless'};

gData.Properties.GroupVariableName = 'ID';

gData.Properties.IndependentVariableName = 'Time';

% Extract variants for RPC_Ag_Tumor from dataset

variants = createVariants(gData, 'RPC_Ag_Tumor', Names='RPC_Ag_Tumor', ...

Model=model, UnitConversion='auto');

variantsForFit = num2cell(variants);

% Create the data dose.

d2 = sbiodose('dose');

d2.TargetName = 'Ab';

d2 = createDoses(gData, 'dose', '', d2);

% Build table of doses.

dosesForFit = num2cell(d2);

% Fitting options

% Define response information.

responseMap = {'RO = RO_exp'};

% Define objects being estimated and their initial estimates.

estimatedInfoObj = estimatedInfo({'kD_Ab'});

estimatedInfoObj(1).Bounds = [1E-1 1E2];

% Define Algorithm options.

options = optimoptions('particleswarm');

options.Display= 'iter';

% Define fit problem.

f = fitproblem('FitFunction', 'sbiofit');

f.Data = gData;

f.Model = model;

f.Estimated = estimatedInfoObj;

f.ErrorModel = 'exponential';

f.ResponseMap = responseMap;

f.Weights = [];

f.Variants = variantsForFit;

f.Doses = dosesForFit;

f.FunctionName = 'particleswarm';

f.Options = options;

f.ProgressPlot = true;

f.UseParallel = false;

f.Pooled = true;

% Run the model with optimum parameters

fitResults = f.fit;

[Mehdi]

Thanks @Jeremy Huard Problem was solved

[Jeremy Huard]

Adding a link to reference page that describes the data format expected by SimBiology as it might be useful to others:

Create Data File with SimBiology Definitions