simFunction and dose Array from an input file

Hi. After fitting my model (fitting was successful using sbiofit) I wrote the code below to generate simfunction to run the model with estimated values and the parmetr in the gData for each case.In gData the doses has been detemined and they are diffent rated doses for each patient applied at diffent time points. I was not able to pass

dosesForFit = createDoses(gData, 'DOSE', 'dose_rate', d3)

to simfunction because

Invalid dosing information. Specify doses as a row vector

and then I used

dosesForFit = reshape(dosesForFit, 1, []); % Ensure dosesForFit is a row vector 

and the problem solved and simFunction was executed but the dose considered the same for all of the cases and was combinatin of all of the doses. How can I solve this problem. I attached input file and the code is below: Thanks

%% Update Parameters with Estimated Values and Running the Model
d3                  = sbiodose('DOSE');
d3.TargetName = 'plasma.A1';
d3.LagParameterName = '';
d3.AmountUnits = 'nanomole';
d3.RateUnits = 'nanomole/day';
d3.TimeUnits=configsetObj.TimeUnits;
dosesForFit = createDoses(gData, 'DOSE', 'dose_rate', d3); 
dosesForFit = reshape(dosesForFit, 1, []); % Ensure dosesForFit is a row vector
doseTable = getTable(dosesForFit);

% getTable(dosesForFit)
% Extract Observables and Experimental Data
for i = 1:length(fitResults.EstimatedParameterNames)
    paramObj = sbioselect(model, 'Name', fitResults.EstimatedParameterNames{i});
    paramObj.Value = fitResults.ParameterEstimates.Estimate(i);
end
% Define parameter names to extract
variantParams = {'BW', 'SLD', 'HSCR', 'FRAPS'};
% Extract IDs in the order they appear in the input file
idsInOrder = exp_data.ID; % Retain original ID order
[uniqueIDs, ~, idOrder] = unique(idsInOrder, 'stable'); % Preserve order from input file
% Initialize paramMatrix with zeros
numIDs = length(uniqueIDs);
paramMatrix = zeros(numIDs, length(variantParams)); % Rows for IDs, columns for parameters
% Populate paramMatrix with parameter values based on the input file order
for i = 1:numIDs
    % Filter data for the current ID
    currentID = uniqueIDs(i);
    idData = exp_data(exp_data.ID == currentID, :);
    % Assign parameter values (assume these are constant per ID)
    paramMatrix(i, :) = [idData.BW(1), idData.SLD(1), idData.HSCR(1), idData.FRAPS(1)];
end
% Observables for the simulation
obs = {'CAb_plasma', 'CADC_plasma', 'CPLun_plasma'}; % Observables to track
% Create SimFunction
sfxn = createSimFunction(model, variantParams, obs, dosesForFit,'UseParallel', true, 'AutoAccelerate', true);
% Simulate Using the Parameter Matrix
% Create dose table from dosesForFit
% Simulate the model using the parameter matrix
simData = sfxn(paramMatrix, [], doseTable, time);

[Arthur Goldsipe]

Can you provide code standalone code that reproduces the problem? As currently shared, the code is not runnable due to several of the variables being undefined. If confidentiality is an issue, would you be able to share the code directly with MathWorks?

[Mehdi]

%% Setup Model
model_input_file='processed_PK_data_DIAG_2.xlsx';
exp_data=readtable(model_input_file);
bystander = 1;
numTumors =7;
setup;
%% Data Import
gData = groupedData(exp_data);
gData.Properties.VariableUnits = {'','','day','day','nanogram/milliliter','nanogram/milliliter','nanogram/milliliter','nanomole','nanomole/day','','kilogram','milliliter','',''};
gData.Properties.GroupVariableName       = 'ID';
gData.Properties.IndependentVariableName = 'TIME';
sbiotrellis(gData,'ID','TIME',{'TAB','CPL','UPL'},'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
variants = createVariants(gData, {'BW','SLD','HSCR','FRAPS'}, Names={'BW','SLD','HSCR','FRAPS'},Model=model, UnitConversion='auto'); % the only thing I need
%% Create the data dose.
d3                  = sbiodose('DOSE');
d3.TargetName = 'plasma.ADC_plasma';
d3.LagParameterName = '';
d3.AmountUnits = 'nanomole';
d3.RateUnits = 'nanomole/day';
d3.TimeUnits=configsetObj.TimeUnits;
dosesForFit = createDoses(gData, 'DOSE', 'dose_rate', d3); 

%% Fitting options
% Define response information.
responseMap = {'CAb_plasma = TAB','CADC_plasma = CPL','CPLun_plasma = UPL'};

% Define objects being estimated and their initial estimates.
estimatedInfoObj = estimatedInfo({'CL','CL_Ab','CL_PL'}); %% for mixed log should be used
% estimatedInfoObj = estimatedInfo({'log(CL)','log(CL_Ab)','log(CL_PL)'}); %% for mixed log should be used

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


rng default


options = optimoptions('particleswarm');
options.Display= 'iter'; 


rng default

% Define fit problem.
f              = fitproblem('FitFunction', 'sbiofit');
% f              = fitproblem('FitFunction', 'sbiofitmixed');

f.Data         = gData;
f.Model        = model;
f.Estimated    = estimatedInfoObj;
f.ErrorModel = 'combined';
f.ResponseMap  = responseMap;
f.Weights      = []; 
f.Variants     = variants;
f.Doses        = dosesForFit;
f.FunctionName = 'particleswarm';
% f.FunctionName = 'nlmefit';

f.Options      = options;
f.ProgressPlot = true;
f.UseParallel  = true;
f.Pooled       = true;
%% Run the model with optimum parameters
fitResults = f.fit;


%% Update Parameters with Estimated Values and Running the Model
dosesForFit = reshape(dosesForFit, 1, []); % Ensure dosesForFit is a row vector
doseTable = getTable(dosesForFit);
% getTable(dosesForFit)
% Extract Observables and Experimental Data
for i = 1:length(fitResults.EstimatedParameterNames)
    paramObj = sbioselect(model, 'Name', fitResults.EstimatedParameterNames{i});
    paramObj.Value = fitResults.ParameterEstimates.Estimate(i);
end
% Define parameter names to extract
variantParams = {'BW', 'SLD', 'HSCR', 'FRAPS'};
% Extract IDs in the order they appear in the input file
idsInOrder = exp_data.ID; % Retain original ID order
[uniqueIDs, ~, idOrder] = unique(idsInOrder, 'stable'); % Preserve order from input file
% Initialize paramMatrix with zeros
numIDs = length(uniqueIDs);
paramMatrix = zeros(numIDs, length(variantParams)); % Rows for IDs, columns for parameters
% Populate paramMatrix with parameter values based on the input file order
for i = 1:numIDs
    % Filter data for the current ID
    currentID = uniqueIDs(i);
    idData = exp_data(exp_data.ID == currentID, :);
    % Assign parameter values (assume these are constant per ID)
    paramMatrix(i, :) = [idData.BW(1), idData.SLD(1), idData.HSCR(1), idData.FRAPS(1)];
end
% Observables for the simulation
obs = {'CAb_plasma', 'CADC_plasma', 'CPLun_plasma'}; % Observables to track
% Create SimFunction
sfxn = createSimFunction(model, variantParams, obs, dosesForFit,'UseParallel', true, 'AutoAccelerate', true);
% Simulate Using the Parameter Matrix
% Create dose table from dosesForFit
% Simulate the model using the parameter matrix
simData = sfxn(paramMatrix, [], doseTable, time);

[Mehdi]

This code doe not contains the setup function for the model, which defines the reactions, species, and other necessary elements. My issue arises after obtaining the estimated parameter values. I want to run the model for all patients and create a SimFunction. However, I cannot pass dosesForFit to the SimFunction without reshaping it into a row vector. When I reshape it to a row vector using reshape, the same dose is applied to all patients, which is incorrect.

I can resolve this issue by running the SimFunction in a loop for each patient, using the following format, but this approach is time-consuming. Is there a more efficient way to solve this?

doseTable = getTable(dosesForFit);

% getTable(dosesForFit)
% Extract Observables and Experimental Data
for i = 1:length(fitResults.EstimatedParameterNames)
    paramObj = sbioselect(model, 'Name', fitResults.EstimatedParameterNames{i});
    paramObj.Value = fitResults.ParameterEstimates.Estimate(i);
end
% Define parameter names to extract
variantParams = {'BW', 'SLD', 'HSCR', 'FRAPS'};
% Extract IDs in the order they appear in the input file
idsInOrder = exp_data.ID; % Retain original ID order
[uniqueIDs, ~, idOrder] = unique(idsInOrder, 'stable'); % Preserve order from input file
% Initialize paramMatrix with zeros
numIDs = length(uniqueIDs);
paramMatrix = zeros(numIDs, length(variantParams)); % Rows for IDs, columns for parameters
% Populate paramMatrix with parameter values based on the input file order
for i = 1:numIDs
    % Filter data for the current ID
    currentID = uniqueIDs(i);
    idData = exp_data(exp_data.ID == currentID, :);
    % Assign parameter values (assume these are constant per ID)
    paramMatrix(i, :) = [idData.BW(1), idData.SLD(1), idData.HSCR(1), idData.FRAPS(1)];
end
% Observables for the simulation
obs = {'CAb_plasma', 'CADC_plasma', 'CPLun_plasma'}; % Observables to track
% Initialize simData to store all results
simData = [];

% Loop through each row of paramMatrix
for i = 1:size(paramMatrix, 1)
    % Extract the parameter set for the current row
    currentParams = paramMatrix(i, :);
    % Create SimFunction
    sfxn = createSimFunction(model, variantParams, obs, dosesForFit(i),'UseParallel', true, 'AutoAccelerate', true);
    % Run the SimFunction for the current parameter set
    currentSimData = sfxn(currentParams, [], doseTable{i}, time);

    % Append the results to simData
    if isempty(simData)
        simData = currentSimData;
    else
        simData = [simData; currentSimData]; % Concatenate results
    end
end

Answer 1

[Arthur Goldsipe]

Ok, I think I finally understand what's going on.

When you create the SimFunction, you only need to specify one "template" dose (for example, d3 or dosesForFit(1)) instead of the entire dose vector. This dose is basically just used to identify the dose target needed for each simulation. In fact, another option would be to just specify the dose target (that is, 'plasma.ADC_plasma') instead of a dose object. You are only required to pass in dose objects (rather than target names) if the dose object has any parameter names specified in its properties (most commonly, lag or duration).

Also, if your goal is to simluate the model using the estimated parameter values, another option would be to use the predict method on the fit results. That is less flexible than creating your own SimFunction, but it's also less work to set up (assuming it meets your needs).

(Sorry if I'm overusing parentheses.)

-Arthur

[Mehdi]

so I canot pass a vector of doses to simFunction? Sould I use sfx in a loop?

[Arthur Goldsipe]

You should use a single dose when creating your SimFunction, but then use a column vector of doses when simulating.

When you create the SimFunction, you only need one dose because you only have one dose target.

When you use the SimFunction, you pass in a column vector of doses. Each row is used in a different simulation.

Let me give one more example. Let's say I want to want to simulating dosing two different species, A and B. I would first create a SimFunction passing in {'A', 'B'} as the dose targets. (Alternatively, I could specify the dose targets by passing in two dose objects with these species as their targets.) The, let's say I want to do 4 simulations with different combinations of dosing. In that case, I'd pass in a 4-by-2 array of dosing information. The first column is the dosing information for species A for each simulation. The second column is the dosing information for species B for each simulation.

Does that clarify things?

[Mehdi]

so how will you pass the simBiology Dose array the first row is for 1 simulatiopn and another one is for the second one. These doses have been read from the input file for differnt patenits which have differnt model parameters. The target of two dose is the same only amount, rate and dosing time is difdfern:

dosesForFit =

SimBiology Dose Array

Index: Name: Type:

1 DOSE (1) repeat

2 DOSE (2) schedule

[Arthur Goldsipe]

After you create the SimFunction, pass in the dose tables as a column vector. That should be how they are returned by createDoses. I think the only change you need to make to your code is to update the createSimFunction line of code to use d3 instead of dosesForFit (although you could also remove the reshape). For completeness, here's an updated excerpt of the code you provided:

d3 = sbiodose('DOSE');
d3.TargetName = 'plasma.ADC_plasma';
d3.LagParameterName = '';
d3.AmountUnits = 'nanomole';
d3.RateUnits = 'nanomole/day';
d3.TimeUnits=configsetObj.TimeUnits;
dosesForFit = createDoses(gData, 'DOSE', 'dose_rate', d3); dosesForFit = createDoses(gData, 'DOSE', 'dose_rate', d3);
sfxn = createSimFunction(model, variantParams, obs, d3.TargetName,'UseParallel', false, 'AutoAccelerate', false);
doseTable = getTable(dosesForFit);
simData = sfxn(paramMatrix, [], doseTable, time);

[Mehdi]

Thanks. Problem solved