Find best input combinations that will maximize the output of a non-linear ODE mathematical model

How to find the best combinations of inputs (3 inputs) that will give the maximum output. Please note that I have a non-linear ODE model with output also as differential eqn. Also, i have SBML file, how to proceed for the same?

[darova]

Look into fmincon

[Darshna Joshi]

@darova@Alan Weiss Below is my matlab code for solving ODE. How to use fmincon for the same. As i am new to it even after looking into the help section of fmincon, couldn't solve the issue.

The problem is to find the optimum(maximum) value of x3 in range of (-8e-4 to 2e-4) by varying kst,x1,x5 and xo)

x5=5 %Input 2 (Input 2 is a state variable and could vary in range of 4 to 15 while performing optimization)

kst=1 %Input 3 (Input 3 is in terms of rate constant, it could vary from 0.1 to 2)

xo=4 %Input 4 (Input 4 is a state variable and could vary in range of 4 to 10)

x1=1e-7 %Input 1 could vary from 1e-9 to 1e-6

%MATLAB CODE%

%Script file

function rest = Scrpt1(t,X)
x2 = X(1); 
x3 = X(2); 

%Parameters

if t<15

x1 = 1e-7; %Input 1 could vary from 1e-9 to 1e-6

else x1 = 0;

end

x5=5 %Input 2 (Input 2 is a state variable and could vary in range of 4 to 15 while performing optimization)

kst=1 %Input 3 (Input 3 is in terms of rate constant, it could vary from 0.1 to 2)

xo=4 %Input 4 (Input 4 is a state variable and could vary in range of 4 to 10)

k1 = 6e7;

km1 = 0.20;

km4 = 0.003;

k3 = 2500.00;

k4 = km4/9;

    km3 = km1;
    
    LAP=1.5
    
    %Differential equations
    dx2dt = km1*x3 + km3*LAP - k1*x1*x2 + km4*x3 - k4*x2;
    
    dx3dt = k1*x1*x2 - km1*(x3+x5+xo) - k3*x3*kst;
    
    rest = [dx2dt; dx3dt];
    
end

% Function file for ODE solution

options = odeset('InitialStep',0.0001,'RelTol',1e-09);
[T,Y]=ode15s(@Scrpt1,[0 60],[9e-13,0],options);
X3= Y(:,2);
plot(T,X3)

How to use fmincon or any other optimization solver for this to solve the mentioned optimization problem of finding maximum value of x3. For which values of x5,kst,xo,x1 we get maximum x3?

My objective function is Maximize(x3).

Answer 1

[Alan Weiss]

Add the parameters to optimize as a vector in your ODE calculation:

function rest = Scrpt1(t,X,params)
x2 = X(1); 
x3 = X(2); 

%Parameters
if t<15
    x1 = params(1); %Input 1 could vary from 1e-9 to 1e-6
else
    x1 = 0;
end

x5 = params(2); %Input 2 (Input 2 is a state variable and could vary in range of 4 to 15 while performing optimization)
kst = params(3); %Input 3 (Input 3 is in terms of rate constant, it could vary from 0.1 to 2)
xo = params(4); %Input 4 (Input 4 is a state variable and could vary in range of 4 to 10)
k1 = 6e7;
km1 = 0.20;
km4 = 0.003;
k3 = 2500.00;
k4 = km4/9;
km3 = km1;
LAP = 1.5;
    
    %Differential equations
    dx2dt = km1*x3 + km3*LAP - k1*x1*x2 + km4*x3 - k4*x2;
    dx3dt = k1*x1*x2 - km1*(x3+x5+xo) - k3*x3*kst;
    rest = [dx2dt; dx3dt];
end

Write your objective function as follows:

function y = objfun(params)
options = odeset('InitialStep',0.0001,'RelTol',1e-09);
[~,Y] = ode15s(@(t,X)Scrpt1(t,X,params),[0 60],[9e-13,0],options);
y = -Y(end,2); % Negative because you want to maximize
% I assume that you want the value of Y at the last time
end

Obtain the solution:

lb = [1e-9,4,0.1,4];
ub = [1e-6,15,2,10];
x0 = [1e-7,5,1,4];
opts = optimoptions("fmincon","PlotFcn","optimplotfval");
[sol,fval] = fmincon(@objfun,x0,[],[],[],[],lb,ub,[],opts)
% The true function value is -fval because you negated y

Alan Weiss

MATLAB mathematical toolbox documentation

[Darshna Joshi]

Thanks a lot sir for such clear explanation. I appreciate your guidance.

[Darshna Joshi]

@Alan Weiss @darovaI used this and now wanted to use the same for SimBiology model. Is it possible to use fmincon for model simulated in SimBiology toolbox?

[Alan Weiss]

Sorry, I don't know anything about SimBiology.

Alan Weiss

MATLAB mathematical toolbox documentation

[Florian Augustin]

Hi Darshna,

you can take a look at sbmimport to create a SimBiology model from your SBML file:

sbModel = sbmlimport("sbmlFileName");

You could add an observable (using addobservable) to the model to compute the value you want to optimize. For example,

% Assuming your model has a component (species, parameter, compartment) Y. 
% funcToComputeQuantityOfInterest is a custom Matlab function you can write. 
% This function takes in the simulated response values of Y and should return
% the scalar value you want to optimize.
addobservable(sbModel, "QOI", "funcToComputeQuantityOfInterest(Y)");

As the next step I suggest to create a SimFunction to simulate your SimBiology model, for example,

% params = names (string vector) of model components in sbModel 
% you want to optimize over.
simFun = createSimFunction(sbModel, params, "QOI", []);

Then you can write your objective function and run the optimization, similar to the workflow Alan described. You can take a look here for how to work with SimFunctions.

function y = objfun(simFun, tStop, params)
[~, y] = simFun(params, tStop);
y = -y{1}(1); % Negative because you want to maximize
% The index (1) is necessary because SimBiology scalar expands the scalar
% observable value to a vector of the same length as the simulation times.
end

% Run the optimization
opts = optimoptions("fmincon", "PlotFcn", "optimplotfval");
tStop = 10; % change to your stop time
params0 = % set your initial values for the optimization; specify the values
          % in the same order as you specified the model components in the
          % params input argument to createSimFunction.
lb = % specify lower bounds
ub = % specify upper bounds
[sol,fval] = fmincon(@(params) objFun(simFun, tStop, params), ...
    params,[],[],[],[],lb,ub,[],opts);

I don't know if it is applicable in your case, but if you are running the optimization for parameter estimation, you could also take a look at SimBiology's sbiofit or sbiofitmixed function.

Best,

Florian

[Darshna Joshi]

Thanks, I will look into that as well.

[Darshna Joshi]

@Alan Weiss I don't know where we can define the range for x3 in which we are supposed to find maximum. As already mentioned in my question The problem is to find the optimum(maximum) value of x3 in range of (-8e-4 to 2e-4) by varying kst,x1,x5 and xo). Although we have got the value of x3 with your code. I feel it is important to mention the range. Please look into that and suggest.

[Alan Weiss]

If the result x3 is inside the range that you want, then there is nothing to do except check it and accept it.

If it is outside the range, then you can include a nonlinear inequality constraint.

function [c,ceq] = nonlcon(x)
% Put the code to calculate x3 here. Then:
ceq = [];
c(2) = x3 - 2e-4; % negative when x3 < 2e-4
c(1) = -8e-4 - x3; % negative when x3 > -8e-4
end

Include @nonlcon in the appropriate place in the fmincon syntax.

Alan Weiss

MATLAB mathematical toolbox documentation

[Darshna Joshi]

Thanks a lot.

[Darshna Joshi]

@Alan Weiss, I tried to solve the problem and the answer was within acceptable limit. However, with different initial search points, the solver is giving different output combinations. Is there any method to choose a good initial search point?

[Alan Weiss]

It is quite difficult to find good start points. I know of no general method. That is one reason that Global Optimization Toolbox offers the MultiStart solver.

Alan Weiss

MATLAB mathematical toolbox documentation

[Darshna Joshi]

@Alan Weiss As per your suggestion, i added nonlcon, but the program is showing error as undefined x3. Although you have included to add the code of x3 here, i am not getting which code am i supposed to add over there in the Function File 1. Is there any error in fmincon regarding the placement of @nonlcon ?

% Function File 1
function [c,ceq] = nonlcon(x)
% Put the code to calculate x3 here. Then:
ceq = [];
c(2) = x3 - 2e-4; % negative when x3 < 2e-4
c(1) = -8e-4 - x3; % negative when x3 > -8e-4
end


% Function File 2
lb = [1e-9,4,0.1,4];
ub = [1e-6,15,2,10];
x0 = [1e-7,5,1,4];
opts = optimoptions("fmincon","PlotFcn","optimplotfval");
[sol,fval] = fmincon(@objfun,x0,[],[],[],[],lb,ub,@nonlcon,opts)
% The true function value is -fval because you negated y

[Alan Weiss]

You need to include the code that accepts the argument x and returns the resulting x3 in the nonlcon function. The nonlcon function has to be able to calculate x3. Then it can impose the constraints.

Alan Weiss

MATLAB mathematical toolbox documentation

[Darshna Joshi]

Okay. Thanks

[Darshna Joshi]

https://in.mathworks.com/matlabcentral/answers/894117-find-best-input-combinations-that-will-maximize-the-output-of-a-non-linear-ode-mathematical-model#answer_781794

@Alan Weiss Will you please guide me that how can i find the global maximum (x3) for the problem mentioned with the use of Genetic Algorithm?

[Alan Weiss]

Sorry, I am not sure that I understand what you are asking. If you are asking how to use ga instead of fmincon as the optimizer, just put ga instead of fmincon. The only difference in syntax is that for ga you give the number of variables in the second argument, nvar instead of an initial point x0. The nvar argument is the length of the x0 vector. Also, ga expects row vectors as arguments, so you might need to change some things, I didn't examine your problem in detail to see whether some of your arguments are matrices or column vectors.

Alan Weiss

MATLAB mathematical toolbox documentation

[Darshna Joshi]

Okay.. Thanks will check it that way @Alan Weiss, I appreciate your timely response always.

Answer 2

[Alan Weiss]

I don't know what an SBML file is. But to optimize an ODE you can look at these examples:

Fit ODE, Problem-Based

Fit an Ordinary Differential Equation (ODE)

Alan Weiss

MATLAB mathematical toolbox documentation

[Darshna Joshi]

@Alan Weiss, i used the code to optimize x3 using fmincon as you stated, i want to apply optimal control to the same, how should i apply optimal control?

[Alan Weiss]

Sorry, I have no idea.

Alan Weiss

MATLAB mathematical toolbox documentation