Borrar filtros
Borrar filtros

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

7 visualizaciones (últimos 30 días)
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);

Respuestas (1)

Jeremy Huard
Jeremy Huard el 6 de Mayo de 2024
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
  4 comentarios
Mehdi
Mehdi el 7 de Mayo de 2024
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
Jeremy Huard el 14 de Mayo de 2024
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;

Iniciar sesión para comentar.

Comunidades de usuarios

Más respuestas en  SimBiology Community

Categorías

Más información sobre Scan Parameter Ranges en Help Center y File Exchange.

Etiquetas

Productos

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by