Global min optimisation for aorta 3-element windkessel model

So, I'm currently trying to optimise a 3-element windkessel model for the aorta using lsqcurvefit (only optimises 1 parameter) and fminsearch (okay, but very dependent on the initial parameters and does not match the pressure or flow exactly or very close). So I am thinking of implementing a global min optimiser, such as the Monte Carlo method but without the large computational cost. Any recommendations?

Answer 1

[William Rose]

@Hussam,

Use fmincon(). I have used it to optimize a model of the circulation with up to 9 parameters, in order to match presure and flow traces recorded from the radial, carotid, and femoral arteries.

A windkessel model, by itself, cannot predict both pressure and flow. It can predict pressure from flow, or flow from pressure, but not both. To get both pressure and flow, you also must have a model of the heart itself, such as a time-varying elastance model. See, for example, this recent article by me and colleagues:

https://onlinelibrary.wiley.com/doi/10.1111/aor.14624

[William Rose]

@Hussam,

Here is an example of data from a simple simulation of a 3-element Windkessel model.

The simulated flow waveform is the positive part of a sine wave, lasting for 40% of the beat duration.

HR=70 bpm. Cardiac output~=5.6 L/min. Mean aortic pressure~=100 mmHg.

Three seconds of simulated data are attached. The columns are time (s), aortic flow (mL/s), aortic pressure (mmHg).

You can fit this with fmincon() to estimate the two resistances and the compliance of the three element model.

data=load('windkesselData.txt');
t=data(:,1); 
q=data(:,2); 
p=data(:,3);
figure; subplot(211), plot(t,q,'-r'); 
grid on; ylabel('Flow (mL/s)'); title('Aortic Flow');
subplot(212), plot(t,p,'-r'); 
grid on; ylabel('Pressure (mmHg)'); title('Aortic Pressure'); xlabel('Time (s)')

OK

[Hussam]

Hi @William Rose, thank you for your feedback. It looks like fmincon could be better than fminsearch, however it is still a local minimum optimizer, correct? Any suggestions for a global min optimizer?

Regarding fmincon(), what parameters/input would you recommend for the function (for the inequalities)?

[Hussam]

I'm trying to use fmincon() as follows, but I get the following error. Kindly advise.

fmincon error nonlcon must be a function

% fmincon inputs
WK_P_ini_fmincon = [5e-10, 1e2, 1e2];
fmincon_lb = [0,0,0];
fmincon_ub = [5e-9, 1e10, 1e10];

Q_peak = 10e-4;  % Peak flow
t_s = 0.28;     % Systolic Time
t_c = 0.84;     % Cycle Time

dt=0.0001;      % (CFD/FSI) Simulation Time Step

%% Importing flow and pressure data

RefQWK = readtable('Ref_P_Q_WK_No-Osc.csv');
RefPWK = readtable('Ref_P_Q_WK_No-Osc.csv');

P_ref=RefPWK.P;                  % multiply by 0.00750062 for mmHg
P0_ref=RefPWK.P(1);              % multiply by 0.00750062 for mmHg

t_ref=RefPWK.t;
t=0:dt:t_ref(end);

P_smooth = smooth(interp1(t_ref, P_ref, t), 50);
P_spline = spline(t,P_smooth);
P_spline_dot = fnder(P_spline);

P=@(t) ppval(P_spline,t);
dPdt=@(t) ppval(P_spline_dot,t);

%Initial P condition for P equation
P0 = P(0);

%Flow equation:
Q_eq = @(t) Q_peak*sin((pi*t)./t_s).^2.*(t>=0).*(t<=t_s);
dQdt_eq = @(t) Q_peak*2*sin(pi*t./t_s).*cos(pi*t./t_s)*pi./t_s.*(t>=0).*(t<=t_s); 

Q = @(t) Q_eq(mod(t,t_c));
dQdt = @(t) dQdt_eq(mod(t,t_c));

%Initial Q condition
Q0 = Q(0);          

%% Curve fitting using the fmincon
 
options = optimoptions('fmincon','Display','iter');

A = [];
b = [];
Aeq = [];
beq = [];

f = @(x) fun_fmin(x, t, P, Q, dQdt, P0);
sol = fmincon(f,WK_P_ini_fmincon,A,b,Aeq,beq,fmincon_lb,fmincon_ub,options);

C_optim_fmincon = sol(1)
R1_optim_fmincon = sol(2)
R2_optim_fmincon = sol(3)

%% Function for fmin optimisation

function f = fun_fmin(x, time, P_in, Q, dQdt, P0) 
    %Renaming the constants
    C = x(1); R1 = x(2); R2 = x(3);
    [t,y] = integration(C,R1,R2,time,Q,dQdt, P0);
    P = y(:,1);
    f = sum((P-P_in(t)).^2);
end

%% Function for ODE Integration using equation

function [t,y] = integration(C,R1,R2,time,Q,dQdt,P0)

  % For Q equation:
  ode_fun = @(t,y)((Q(t)/C)*(1+R1/R2) + R1*dQdt(t) - (y(1)-0)/(C*R2));
  
  [t,y] = ode45(ode_fun, time, P0);
  
end

[William Rose]

@Hussam,

If you desire, you can explore the tools in the global optimization toolbox. Simulated annealing, particle swarm, and multiple starting point algoritms may be worth considering. The only method I have used is multiple starting points. There is now a toolbox item, multistart (here and here), to assist in this process. I started using multiple starting points before there was a global optimization toolbox, so I just wrote my own implementation (see below), and it seems to work fine. I have used it successfully to estimate parameters for cardiovascular models with up to nine parameters (resistors, capacitors, inertances, wave reflection).

You may decide not to use Matlab's multistart. There is a bit of a learning curve. It's not hard to do it yourself. Just call fmincon() with various initial guesses for the 3 parameters, and pick whichever solution gives the smallest error. For a three parameter model, the search space is a 3D cube, bounded by the lowest and highest reasonable values for each parameter. You could start at 9 points: the center of the cube, and 10% or 20% in from each of the 8 corners. Or you could start at 27 points: 10%, 50%, and 90% along each of the three axes.

For fmincon(), you don't need to use the equality or inequality constraints. You can just use [ ] for those arguments, when calling fmincon(), a shown:

x = fmincon(fun,x0,[],[],[],[],lb,ub);

Good luck.

[William Rose]

@Hussam,

YOu need a [] before the options, when passing parameters to fmincon.

In other words, you need

sol = fmincon(f,WK_P_ini_fmincon,A,b,Aeq,beq,fmincon_lb,fmincon_ub,[],options);

I don;t know if that will make everything work, but it should solve that one error.

Your approach is complicated. I see that you have created functions P(t), Q(t), dPdt(t), etc. I'm not sure why you don;t use a vector of pressures and a vector of flows. I also see that you interpolate the simulation results to a time base with dt=0.0001. What are the units for time, and pressure, and flow? If time is in seconds, then I don;t see why you need such a high sampling rate. I think 200 Hz would suffice. This would reduce the computational load.

[Hussam]

Thank you for the feedback. I got it to work with MultiStart and with the SQP algortihm and the interior-point algorith. However, whenever I try to add 15 or more points, the optimizer gets stuck after a few points. Any advice on this?

Also, I am setting the LB as zero, would that be a problem? I see that while it iterates that some objective functions return NaN (maybe because it is trying C = 0?, not sure), but it deals with that anyway by going for another point.

[William Rose]

@Hussam,

I have not used MultiStart so I cannot advise you about it.

I am not sure what you mean when you say "whenever I try to add 15 or more points...". Do you mean that the pressure and flow signals must be 14 points long or shorter? That would be strange. I have fitted signals that are hundreds or thousands of points.

When you say LB=0, I assume LB refers to lower bound, and I assume you mean LB=[0,0,0], because your model has three parameters. I would not use LB=[0,0,0], due to the possibilty of NaNs. I would specify very small positive values for the lower bound vector. For the resistances, I would specify a lower bound that is 1% of total peripheral resistance (TPR=mean pressure/mean flow). In your model data, the TPR is about 95 mmHg/(80 mL/s) ~= 1.2 mmHg-s/mL. (Convert to whatever units you are using.)

For C, notice that the time constant of decay of arterial pressure, when there is no flow, is on the order of 1 second. Therefore a very rough initial guess for C is C*TPR=1 s, i.e. C=(1 s)/TPR ~= 0.8 mL/mmHg. Choose a LB for C that is 1% of 0.8 mL/mmHg.

I plotted your P and F signals. At the end of diastole, the pressure is about 80 mmHg and is not changing at all. The fact that the "final" pressure asymptote is about 80 mmHg is not consistent with a standard 3 element windkessel model. In a standard 3 element windkessel, the pressure will decay to zero if you wait a long time. Therefore I do not expect a 3 element windkessel to fit your data well.

If you have more questions, please email me securely by clicking on the envelope icon in the pop-up window that appears when you click on the "WR" next to my posts. I may not answer for about two weeks, since I am very busy starting tomorrow.

[Hussam]

Hi @William Rose, thank you very much for your feedback. With number of points, I meant the number of points for multistart. Sometimes fmincon() gets stuck with certain points.

Regarding C I did put down 1e-10 m^3/Pa (where your suggested calculcation is 6e-9, not very far off).

Regarding my pressure decay, would you recommend me to use a higher order model (4 element for example), or should I try to have a pressure waveform that decays to zero. What do you think? Thank you!

[William Rose]

@Hussam,

You're welcome.

"Would you recommend me to use a higher order model (4 element for example), or should I try to have a pressure waveform that decays to zero. What do you think?"

To answer your question about how to proceed (for example, use a 4-4lement Windkessel), I would like to understand where the P and F data came from*, and what your goals are. Please email me securely by clicking on the envelope icon in the pop-up window that appears when you click on the "WR" next to my posts. I may not answer for a day or two, due to other obligations.

* I see "CFD/FSI" in a comment in your code. Did your data come from a CFD simulation involving fluid-structure interaction?