System identification Windkessel - pulse wave

Hello, I would like to know if there is any way how to identify parameters if I have specified input to model and desired output. To be more specific I am modelling Windkessel model. It´s electrical model from two resistors one inductor and one capacitor. Wanted output is pulse wave.

Model:

And this is desired output of model (voltage at resistor R):

Can anyone help me please? Thank you

[Rodrigo Vialle]

Hello, could you share your latest model?

Answer 1

[Star Strider]

If you provide the input and output as a ‘.mat’ file, and your model, it should be possible to derive the component values. (I can derive a circuit model for it if you want.)

You aren’t doing system identification since you already have a model for your system. You’re doing parameter estimation of your model.

EDIT —

The circuit model I derived for your model (I let the Symbolic Math Toolbox do the algebra) is:

syms C L R r Pi Po P2 s t 
Z1 = 1/r;
Z2 = s*C + 1/R;
Z3 = 1/(s*L);
Zt = Z1 + Z2 + Z3;

I = Pi/Zt;
Po = Z2*I;
Out = ilaplace(partfrac(Po, s));
Out = simplify(Out, 'Steps',20);
Out_fcn = matlabFunction(Out, 'Vars',{t,Pi,r,C,R,L})

Out_fcn = @(t,Pi,r,C,R,L) Pi.*dirac(t)-(Pi.*exp((t.*(R+r).*(-1.0./2.0))./(C.*R.*r)).*(cosh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r))-1.0./sqrt(L).*sinh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r)).*(L.*r+L.*R-C.*R.*r.^2.*2.0).*1.0./sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0)))./(C.*r);

The ‘Z’ values are the three lumped impedances (and ‘Zt’, the total impedance), where the output is taken across ‘Z2’. The idea is straightforward: the current ‘I’ (or blood flow here) through the series of impedances is simply ‘Pi/Zt’. It is then straightforward to calculate the voltage (or pressure) drop across ‘Z2’ as ‘Z2*I’.

In order to fit your function and estimate the parameters, we need three data vectors: ‘t’, ‘Pi’, and ‘Po’, the time, input pressure with respect to time, and output pressure with respect to time. The function fits the output pressure it estimates with specific parameters to the ‘Po’ vector.

The procedure for estimating the parameters would then be:

% % % MAPPING: b(1) = r,  b(2) = C,  b(3) = R,  b(4) = L
Out_fcn = @(t,Pi,r,C,R,L) Pi.*dirac(t)-(Pi.*exp((t.*(R+r).*(-1.0./2.0))./(C.*R.*r)).*(cosh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r))-1.0./sqrt(L).*sinh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r)).*(L.*r+L.*R-C.*R.*r.^2.*2.0).*1.0./sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0)))./(C.*r);

t = t(:);                                                           % Column Vector
Pi = Pi(:);                                                         % Column Vector
Po = Po(:);                                                         % Column Vector
b0 = ['Initial Parameter Estimates For r,C,R,L — IN THAT ORDER — As A Vector'];
nrmrsd = @(b) norm(Po - Out_fcn(t,Pi,b(1),b(2),b(3),b(4)));         % Norm Of Residuals (Errors)
bvals = fminsearch(nrmrsd, b0);

The ‘bvals’ vector will be the estimated parameters.

NOTE — This is UNTESTED CODE since I do not have your data to test it with.

[David Hatle]

Hello, thank you very much I really appreciate your help. I provide input and output data and simulink model. If you help me with that it would be almost lifesaving! There is a little change in model arrangement - inductor is connected parallely to resistor r (showed in WK.slx) sorry for that. Thank you again!

[Star Strider]

My pleasure!

I’ll put ‘r’ in parallel with ‘L’ and let the Symbolic Math Toolbox re-derive the model. I’ll post it here.

I can’t guarantee success with the parameter estimation, because actual physiological systems are notoriously difficult to identify.

I also need a time vector, so if you didn’t provide one (I’ve not loaded the 'mat files yet), I need to know the sampling interval (time in seconds between samples, assuming a constant sampling interval). I assume the lengths are the same for both ‘input’ and ‘output’, so the duration should be straightforward once I have the sampling interval.

I have done modeling with various physiological systems including Windkessel models (I’m a Board Certified Internal Medicine Physician and M.S. Biomedical Engineer). However I need for you to provide approximate initial values for the parameters ‘r,C,R,L’.

I will re-derive the model, and wait for you to provide the approximate initial parameter estimates. Without those, there could be any number of sets of parameter estimates that could fit, all but one of which would be completely meaningless to you, or that could result in a complete failure to fit your data. That’s simply the nature of nonlinear parameter estimation. So I very definitely need to know what the initial estimates are.

EDIT —

The re-derived model is:

Out_fcn = @(t,Pi,r,C,R,L) (Pi.*exp((t.*(R+r).*(-1.0./2.0))./(C.*R.*r)).*(cosh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r))-1.0./sqrt(L).*sinh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r)).*(L.*r+L.*R-C.*R.*r.^2.*2.0).*1.0./sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0)))./(C.*r);

All I now need are your initial parameter estimates and we’re cleared to go!

[David Hatle]

There are two time vectors.. both durations are 1 second but number of samples (lengths) in both vectors are different. Should I resample them?

To be honest I don´t know approximate inital values. I made a script that compared reference pulse wave with wave created by windkessel. Comparison according to correlation worked quite well but I also had problem with initial values and algorithmic complexity was nearly polynomial so when trying combination of multiple parameters in for cycle it was a dead end.

Is there any way how to do it without initial estimation? I did some experiments and it was a little similar to pulse wave so let´s give it a try: r = 50 Ω; C = 4e-4 F; L = 1.5e3 H; R = 1050 Ω;

I am also student of biomedical engineering :)

[Star Strider]

The time vectors both have to be the same. I can use the resample function (preferably) or interp1 (last resort for signal processing) to put them on the same time reference. They won’t match (there are phase and other considerations), but they have to have the same time base and length or MATLAB will throw a deservedly unrecoverable error.

If you know that the sampling interval is 1 sec, then I would truncate the longer sequence to the length of the shorter sequence. I will do so in my code.

Please note that 1.0s is a very long sampling interval in circulatory dynamics, since the pulse interval of humans is about 0.8s or less. I would choose a sampling interval of 0.5s, and preferably 0.1s, for circulatory dynamics studies.

I know it consumes memory and disk space. It also results in publishable data. Err on the side of publishable data. That’s why whatever Deity you believe in (or don’t) created DVD disks and 64 GB micro SD cards.

The initial parameter estimates are important, since I don’t know the system you’re studying. Values differ with respect to age and species. Yours will provide an acceptable initial estimate for the system you’re studying.

[David Hatle]

I have specified system I am studying but it says nothing about initial values. I have pulse wave measured on radial artery (output) and data that shloud approximate ejection of left ventricle - positive sine wave (input). Input should simulate ejection fraction which takes app 0.2 s. But I have no idea what resistance, compliance or inertia of this part of arterial tree might be. I would have one more question if you don´t ming me asking... how did you obtain Out_ fcn? I know it´s based on those parameters but how to get if from provided time series?

[Star Strider]

I may use the Global Optimization Toolbox functions to find the best parameters. That could take a while.

I derived the ‘Out_fcn’ (in R2016b) from the information in your revised circuit diagram descriptions:

syms C L R r Pi Po P2 s t 
Z1 = 1/(1/r + 1/(s*L));                                     % Parallel ‘rL’ Network
Z2 = 1/(s*C + 1/R);                                         % Parallel ‘RC’ Network
Zt = Z1 + Z2;                                               % Total Impedance

I = Pi/Zt;                                                  % Network Current
Po = Z2*I;                                                  % Output Pressure
Out = ilaplace(partfrac(Po, s));                            % Inverse Laplace Transform
Out = simplify(Out, 'Steps',20);
Out_fcn = matlabFunction(Out, 'Vars',{t,Pi,r,C,R,L})

The rest of my code remains the same.

The radial artery information helps. I may go online with PubMed to see if anyone else has identified relevant parameters.

[Star Strider]

I opted to interpolate the longer sequence to the shorter 201-element sequence. That way, I avoid creating data that aren’t actually present in the shorter sequence.

That aside, we have a problem with the data. Specifically, when I load your data, I notice some significant differences in the magnitudes of the source and output pressures:

inp = load('David Hatle input.mat');
otp = load('David Hatle output.mat');

Po = interp1(otp.t_output, otp.x_output, inp.time);

q1 = [inp.src(1:10)'  Po(1:10)']

q1 =
           0  8.9981e+03
  2.1565e-02  8.9660e+03
  4.2998e-02  8.9357e+03
  6.4166e-02  8.9064e+03
  8.4938e-02  8.8760e+03
  1.0519e-01  8.8446e+03
  1.2479e-01  8.8144e+03
  1.4362e-01  8.7866e+03
  1.6157e-01  8.7596e+03
  1.7852e-01  8.7314e+03

The input pressures are in the first column, and the output pressures are in the second column. The network (that I believe I have derived correctly) is purely passive. It does not amplify.

From a purely technical perspective, another problem is that the source data begins at zero. That simply does not occur in normal physiology. The systemic arterial blood pressure limits are characteristically between about 80 mm Hg (Torr) or 10.5 kPa to 120 mm Hg or 16 kPa in a normal, healthy human. (The diastolic can be as low as 40 mm Hg or 5.3 kPa with aortic valve regurgitation, since the pressures are significantly lower in the pulmonary circulation, and left ventricular pressure has to drop below the pulmonary vein pressure in order to fill.) Aortic root pressure values in the range you have them are incompatible with life.

What are the units of the output pressures? Pascals?

I will rely on you to make the necessary corrections and scale the pressures appropriately. It is necessary to get the pressure units correct in order for the identified component values to have any meaning. Then we can try this again.

[David Hatle]

Yes right sorry that´s my bad. I scaled data and name it with respect to notation you use. Pressures are in mmHg. I attach modified data. Thanks again for your effort and patience.

[Star Strider]

My code now estimates the parameters without returning a NaN value because of the initial zero value for the pressure. It still doesn’t give anything close to an appropriate fit. I even used patternsearch to see if it could find a better parameter set, but the result was essentially the same as with fminsearch.

I also tried shifting the aortic pressure waveform to approximately coincide with the radial artery waveform (see the Wikipedia article on the Cardiac cycle), again without success.

I used this code to do the shift:

idx = find(t <= 0.1, 1, 'last');
Pi = circshift(Pi, idx+1);

The problem may be your model. I went online to PubMed (using this search string) and found 18 articles that looked at Windkessel models in radial artery circulatory dynamics. I found one slightly informative (and free) article Noninvasive Pulse Wave Analysis for the Early Detection of Vascular Disease. Hypertension 1995;26:503-508. It’s not recent, but you can use its PMID: 7649589, as well as the search string, to search for others. (Most of my Windkessel experience involved ‘designing’ various pulse waveforms to be programmed into a pulsatile perfusion pump for cardiopulmonary bypass during my brief but spectacular foray into private industry long, long ago in a galaxy far, far away. My model belonged to my employer, so I wasn’t able to take it with me. It wasn’t published outside of our successful FDA Class III PMA application for the pump, and would probably be redacted anyway should anyone request it.) None of my Windkessel experience is recent.

I’ll be glad to work with you on this as my time permits. It’s certainly an interesting problem.

[David Hatle]

I will definitely do more research about this topic.

The pleasure is mine! I would be greatful to work with so skillful programmer like you. I am sure you are much more experienced than me. On the other hand I have a great motivation beacuse my master degree depends on it :). Would you be willing to give me some contact on you (maybe an email) that would make our conversation easier? Thank you for reply

[Star Strider]

My pleasure.

I generally do not give help other than here on MATLAB Answers. (The exceptions are only to resolve minor ambiguities that do not affect the substance of the Answers threads.) To do otherwise destroys the evolution of problem solutions, and creates discontinuities that others searching for help on the same (or similar) topics would need to solve their own problems.

My M.S. involved attempting to model the regulatory dynamics of carbohydrate metabolism. (I am Certified by the American Board of Internal Medicine, did a fellowship in Endocrinology.) That problem turned out to be much more complicated than I ever imagined, and is a problem that remains unsolved because of the complexity of the physiological and biochemical processes involved. I got ‘close enough’ to satisfy my committee. I am convinced that they awarded me my degree at least in part because they no longer wanted to hear my screams of frustration.

Physiological regulatory processes (one of which you are modeling) evolved over millions of years to be adaptive to myriad external disturbances in order to keep us alive, providing endless M.S. and Ph.D. research opportunities in the process (and endless frustration for clinicians like me when things go awry), since there is more to them than will ever be dreamt of in any of our philosophies.

That is the reason I find your research personally interesting. There is much to be discovered, and it is quite possible that you could find something new.

Physiological processes are in general nonlinear, nonstationary, and nontrivial. Enjoy the challenge!

[David Hatle]

Undestood. I am trying now to get data into Simulink Design Optimization toolobox which should do exactly what I want. But I am having some trouble with inport and output data.

[Star Strider]

I have the Toolbox, but I have no experience with it.

Simply reading the documentation, note that it first requires an accurate model, and that the model already be present as a Simulink model, so it seems you have to do that first. It then will allow you import the data and to tweak the parameters. You quite probably need to explore alternate models, preferably one that comes close to reproducing your radial pulse waveform.

Some Windkessel models that I have seen incorporate a rectifier or gated rectifier. It would likely be best for you first to explore the literature to see what others have done. That is where I always begin.

[David Hatle]

I made the Simulink Design Optimization toolbox work but the result is not satisfying. Due to your experience with windkessels - do you know if it is possible to model curve with simple 4-element windkessel (my model with 4 components) that has dicrotic notch? Or it is impossible to have such specific curve generated by this simple model. Thank you

[Star Strider]

My pleasure.

As I remember, modeling the dicrotic notch requires a rectifier of some sort in order to simulate the closing of the aortic valve. I do not recall the details at present. (A rectifier with a bias can be a challenge to model. I do not remember now how I did it before.)

You can get some information from a 4-component model, as did the authors of the paper I referenced. It depends on what you want from the model.

[David Hatle]

With some model adjustments I was able to find parameters for all components. And these are the normalised results:

Now I will focus on modelling dicrotic notch.

[Star Strider]

The ‘model adjustments’ must make your new model significantly different from the model you first posted?

Your posted model could not produce this latest result.

[David Hatle]

Yes... I used current source instead of voltage source so the amplitude was much smaller. Then I extended model to circuit of 8 components - "divided" into 3 parts (chest/abdomen, arm and leg). Then I let Simulink Design Optimization do the rest. Of source initial estimation of all parameters was necessary as we discued earlier. This waveform is taken from resistor representing arm. Desing optimization toolbox is really migthy tool for parameter estimation.

[Star Strider]

Reconfiguring the model, adding 4 parameters, certainly are as well.

I assume you no longer need my help.

[David Hatle]

Well the question has been answered so I will accept your answer. Thank you. My work is definitely not finished at this point... there is still much work to do. Do you think if is it possible to contact you in case of some hard to solve problem concerning this matter?

[Jan]

@David: Please press the "accept" button to accept the answer.

[Star Strider]

@David — If you still need my help, I will provide it as much as I am able. I prefer to respond on MATLAB Answers to provide continuity. I simply felt ‘out of the loop’ with your new model (that you have not shared here).

@Jan — Thank you!

[David Hatle]

Dear Star Strider I would have another question... I have had succes with current source but now I am trying to simulate state when I have pressure input to controlled voltage source and I measure voltage at resistor R (scheme at the top of this topic). I can´t figure out how to convert units from fluid domain to electric domain. I found really good article were the fluid units are described. For example capacitance is in cmˆ3/mmHg or resistance in mmHg/cmˆ3/ sec. Is there any way how to calculate them to farads and ohms so I can use simulink? Thank you

[Star Strider]

I am not aware of any way to convert units of liquid pressure, flow, and resistance to electrical units. (The exception being calibrated transducers where a transducer would be calibrated in terms of volt/torr or some such.) They are completely different physical properties. You would probably have to go back to the original definitions of the relevant units to see if you could determine direct relationships between them. (Blood has a critical Reynolds number of 2000 if I remember correctly, and while it generally stays below that value, blood is non-Newtonian and has nonlinear characteristics.)

I have never seen any conversions, but they may exist. I am not certain how to search the literature or the Internet for such information, although a physiology reference text could have them. A brief search found Units and Conversion Factors and List of equations in fluid mechanics. A more extensive search could discover what you want.

Answer 2

[John BG]

Hi David

1.

I have simulated the circuit of your question in ADS and I get the following transient

2.

tuning C1 and L1 both min values

3.

C2 and L1 both max values

4.

R2 value is quite irrelevant, the one that counts is R1/R2 it has to be kept small

5.

In MATLAB you can either build a model from ground zero, or could start with Simulink. Do you have Simscape?

So, David

if you find this answer useful would you please be so kind to mark my answer as Accepted Answer?

To any other reader, please if you find this answer of any help solving your question,

please click on the thumbs-up vote link,

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

[David Hatle]

Thank you for your answer but to be honest I am not wiser after reading your answer. Yes I have simscape. I found out that this task can be acomplished by using genetic algorithm... are you familiar with that? Thank you. FYI I attach input, output and simulink model. Thank you for your help

[John BG]

1.

I use MATLAB as often as possible, but for the task you asked to solve and for the time I allocated to your question, ADS save time, even if using Simscape.

2.

So, Simscapte it is, would you please attach to your question your Simscape circuit?

3.

Also, do you want to generate an ECG signal, or the wavelet transform of an ECG?

4.

Would the following in anyway help?

the Illustrated wavelet transform, by Paul Addison, CRC press, page 297

https://books.google.co.uk/books?id=wBoNDgAAQBAJ&pg=PA275&lpg=PA275&dq=heartbeat+time+template+standard&source=bl&ots=ybcnc8_Vk8&sig=V5VOO8zsTDNdadOKNzrbRZZdUfI&hl=en&sa=X&ved=0ahUKEwiSidHtsb7SAhXnI8AKHXkRDd4Q6AEIVjAJ#v=onepage&q=heartbeat%20time%20template%20standard&f=false

[David Hatle]

Needed data including Simscape model are under this comment. ECG is not specifically required. I only need to get parameters of circuit model based on input and output data.