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

Question

Mehdi 2024-4-25

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2111806-sbiofit-ability-to-estimate-the-parameter-for-3-different-cases-at-the-same-time

评论： Jeremy Huard 2024-7-3

采纳的回答： Jeremy Huard

exp.xlsx

在 MATLAB Online 中打开

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);

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jeremy Huard 2024-5-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2111806-sbiofit-ability-to-estimate-the-parameter-for-3-different-cases-at-the-same-time#answer_1453121

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Mehdi 2024-5-7

在 MATLAB Online 中打开

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 2024-5-14

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 2024-6-20

Thanks @Jeremy Huard Problem was solved

Jeremy Huard 2024-7-3

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

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

社区

另请参阅

类别

标签

产品

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

社区

另请参阅

类别

标签

产品

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论