Shifting an array to create a time delay within a function.

I am doing a project that involves modelling the cardiovascular system. I have a function that generates the system of differential equations that make up the model, and this function is fed into the ode23t solver. In the model, the ventricles and atria are controlled by separate 'activation functions'. These are time functions that control when the chambers contract. An important feature of the model is that the atria contract before the ventricles. Therefore, I want the atrial activation function to depend on one time array and the ventricular activation function to depend on another, shfted time array. This is how I am trying to implement it in the function:

function dx1=CVS(t1,x1)

% set-up (excluded here for brevity)

%set timing variables
HR = 70; % heart rate (bpm)
HR_s = HR/60; % heart rate (bps)
T = 1/HR_s; % period of cardiac cycle (s)
T_max_a = 0.125; % time to end systole for atria (s)
T_max_v = 0.2; % time to end systole for ventricles (s)

tm1 = mod(t1,T);
tf1 = 5;

tau_a = 0.025; % time constant of relaxation for atria (s)
tau_v = 0.03; % time constant of relaxation for ventricles (s)

%calculate AV delay in terms of increments in the t1 array

AV_delay = 0.16; %AV delay (s)
interval_length_s = tf1; %length of interval (s)
t_per_s = length(t1)/interval_length_s; %increments in t1 per second
shift_t = ceil(AV_delay*t_per_s); %AV delay in terms of t1 increments

%delay tm
n = shift_t;
tm2 = circshift(tm1,n);
if n>0
    tm2(1:n) = 0;
else
    tm2(end+n+1:end) = 0;
end

%activation function for atria
if tm1 < 1.5*T_max_a
    a1 = 0.5*(1 + sin((pi*tm1/T_max_a) - (pi/2)));
else
    a1 = 0.5*exp(-(tm1 - 1.5*T_max_a)/tau_a);
end

%activation function for ventricles

if tm2 < 1.5*T_max_v
    b1 = 0.5*(1 + sin((pi*tm2/T_max_v) - (pi/2)));
else
    b1 = 0.5*exp(-(tm2 - 1.5*T_max_v)/tau_v);
end
 
end

And this is where the function is called:

t01 = 0;
tf1 = 5;
interval1 = t01:0.01:tf1;

%initial conditions
x01 = [75.6 69 65 118.6 72.7 61.2 408 79.8];

[t1,x1]=ode23t(@CVS, interval1, x01);

I have tested this outside the function and it works as expected. However, the activation functions inside the function do not behave in the correct way.

Is there a way to implement this kind of time shift inside the function?

Any help is appreciated.

Answer 1

[Torsten]

Use dde23 instead of ode23t.

[Hannah Goldblatt]

Thank you for your response. I'm not sure if that would work. The ODEs are not time-shifted versions of each other. It's more like I have a series of equations, some of which contain time functions. This is an example of the type of system I am modelling:

The working code for the model is included below. I would like to shift this function 16ms forward:

tm1 = mod(t1,tc1);

%impulse fxn
if tm1 < 1.5*T_es
    a1 = 0.5*(1 + sin((pi*tm1/T_es) - (pi/2)));
else
    a1 = 0.5*exp(-(tm1 - 1.5*T_es)/tau);
end

I tried to do it like this:

%delay tm
n = shift_t;
tm2 = circshift(tm1,n);
if n>0
    tm2(1:n) = 0;
else
    tm2(end+n+1:end) = 0;
end

%impulse fxn
if tm2 < 1.5*T_es
    a1 = 0.5*(1 + sin((pi*tm2/T_es) - (pi/2)));
else
    a1 = 0.5*exp(-(tm2 - 1.5*T_es)/tau);
end

But replacing tm2 with the shifted version doesn't seem to work.

Here is the working code for the model:

clc;
clear;

V_0 = 0; % unstressed volume (ml)
A = 0.35; % scaling factor for EDPVR (mmHg)
B_lv = 0.033; % exponent for LV EDPVR (1/ml)
B_rv = 0.023; % exponent for RV EDPVR (1/ml)
E_es_lv = 3; % LV end-systolic elastance (mmHg/ml)
E_es_rv = 0.7; % RV end-systolic elastance (mmHg/ml)

%Circulation

%Pulmonary
R_ap = 40/1333.33; % arterial resistance (mmHg.s/ml)
R_cp = 27/1333.33; % characteristic resistance (mmHg.s/ml)
R_vp = 20/1333.33; % venous resistance (mmHg.s/ml)

C_ap = 13; % arterial capacitance (ml/mmHg)
C_vp = 8; % venous capacitance (ml/mmHg)

%Systemic
R_as = 1200/1333.33; % arterial resistance (mmHg.s/ml)
R_cs = 40/1333.33; % characteristic resistance (mmHg.s/ml)
R_vs = 20/1333.33; % venous resistance (mmHg.s/ml)

C_as = 1.32; % arterial capacitance (ml/mmHg)
C_vs = 70; % venous capacitance (ml/mmHg)

t01 = 0;
tf1 = 4;

%interval in which CVS behaviour is evaluated

interval1 = [t01 tf1];

%initial conditions

x01 = [102.5 91.5 83.1 245.8 145.7 81.2];

%specify ode23t as resolution method for ODEs

[t1,x1]=ode23t(@CVS,interval1,x01);

%from time matrix (t) and state variables (xi, i = 1,...,6), create the
%non-state variables (pressures), which cannot be plotted directly by
%ode23t. 

HR = 75; % heart rate (bpm)
HR_ms = HR/60; % heart rate (bps)
tc1 = 1/HR_ms; % period of cardiac cycle (s)
T_es = 175/1000; % time to end systole (s)
T_max = 0.2;
tau = 25/1000; % time constant of relaxation (s)

tm1 = mod(t1,tc1);

a1=zeros(length(t1),1);

for k = 1:length(t1)
    if tm1(k) < 1.5*T_es
        a1(k) = 0.5*(1 + sin((pi*tm1(k)/T_es) - (pi/2)));
    else
        a1(k) = 0.5*exp(-(tm1(k) - 1.5*T_es)/tau);
    end
end

P_es_lv = E_es_lv*(x1(:,1) - V_0); % LV ESPVR
P_ed_lv = A*(exp(B_lv*(x1(:,1) - V_0)) - 1); % LV EDPVR
P_es_rv = E_es_rv*(x1(:,2) - V_0); % RV ESPVR
P_ed_rv = A*(exp(B_rv*(x1(:,2) - V_0)) - 1); %RV EDPVR

P_lv = (P_es_lv - P_ed_lv).*a1 + P_ed_lv; %left ventricular pressure
P_rv = (P_es_rv - P_ed_rv).*a1 + P_ed_rv; %right ventricular pressure

%systemic and pulmonary circulation

P_C_as = x1(:,3)/C_as; % systemic arteries
P_C_vs = x1(:,4)/C_vs; % systemic veins
P_C_ap = x1(:,5)/C_ap; % pulmonary arteries
P_C_vp = x1(:,6)/C_vp; % pulmonary veins

%Kirchoff's laws for diodes. Simulation of valves.

for k = 1:length(a1)

    if P_lv(k) < P_C_vp(k)
        alpha_lv(k) = 1;
    else
        alpha_lv(k) = 0;
    end

    if P_rv(k) < P_C_vs(k)
        alpha_rv(k) = 1;
    else
        alpha_rv(k) = 0;
    end

    if P_lv(k) > P_C_as(k)
        beta_lv(k) = 1;
    else
        beta_lv(k) = 0;
    end

    if P_rv(k) > P_C_ap(k)
        beta_rv(k) = 1;
    else
        beta_rv(k) = 0;
    end

end

figure
tiledlayout('flow')
nexttile
plot(t1,x1(:,1)); title('LV Stressed Volume'); ylabel('Volume (ml)'); xlabel('Time (s)');
nexttile
plot(t1,x1(:,2)); title('RV Stressed Volume'); ylabel('Volume (ml)'); xlabel('Time (s)');
nexttile
plot(t1,x1(:,3)); title('Systemic Arterial Stressed Volume'); ylabel('Volume (ml)'); xlabel('Time (s)');
nexttile
plot(t1,x1(:,4)); title('Systemic Venous Stressed Volume'); ylabel('Volume (ml)'); xlabel('Time (s)');
nexttile
plot(t1,x1(:,5)); title('Pulmonary Arterial Stressed Volume'); ylabel('Volume (ml)'); xlabel('Time (s)');
nexttile
plot(t1,x1(:,6)); title('Pulmonary Venous Stressed Volume'); ylabel('Volume (ml)'); xlabel('Time (s)');

function dx1=CVS(t1,x1)

V_0 = 0; % unstressed volume (ml)
A = 0.35; % scaling factor for EDPVR (mmHg)
B_lv = 0.033; % exponent for LV EDPVR (1/ml)
B_rv = 0.023; % exponent for RV EDPVR (1/ml)
E_es_lv = 3; % LV end-systolic elastance (mmHg/ml)
E_es_rv = 0.7; % RV end-systolic elastance (mmHg/ml)

%Circulation

%Pulmonary
R_ap = 40/1333.33; % arterial resistance (mmHg.s/ml)
R_cp = 27/1333.33; % characteristic resistance (mmHg.s/ml)
R_vp = 20/1333.33; % venous resistance (mmHg.s/ml)

C_ap = 13; % arterial capacitance (ml/mmHg)
C_vp = 8; % venous capacitance (ml/mmHg)

%Systemic
R_as = 1200/1333.33; % arterial resistance (mmHg.s/ml)
R_cs = 40/1333.33; % characteristic resistance (mmHg.s/ml)
R_vs = 20/1333.33; % venous resistance (mmHg.s/ml)

C_as = 1.32; % arterial capacitance (ml/mmHg)
C_vs = 70; % venous capacitance (ml/mmHg)

HR = 75; % heart rate (bpm)
HR_ms = HR/60; % heart rate (bps)
tc1 = 1/HR_ms; % period of cardiac cycle (s)
T_es = 175/1000; % time to end systole (s)
T_max = 0.2;
tau = 25/1000; % time constant of relaxation (s)

tm1 = mod(t1,tc1);

if tm1 < 1.5*T_es
    a1 = 0.5*(1 + sin((pi*tm1/T_es) - (pi/2)));
else
    a1 = 0.5*exp(-(tm1 - 1.5*T_es)/tau);
end

P_es_lv = E_es_lv*(x1(1) - V_0); % LV ESPVR
P_ed_lv = A*(exp(B_lv*(x1(1) - V_0)) - 1); % LV EDPVR
P_es_rv = E_es_rv*(x1(2) - V_0); % RV ESPVR
P_ed_rv = A*(exp(B_rv*(x1(2) - V_0)) - 1); %RV EDPVR

P_lv = (P_es_lv - P_ed_lv)*a1 + P_ed_lv; %left ventricular pressure
P_rv = (P_es_rv - P_ed_rv)*a1 + P_ed_rv; %right ventricular pressure

%systemic and pulmonary circulation

P_C_as = x1(3)/C_as; % systemic arteries
P_C_vs = x1(4)/C_vs; % systemic veins
P_C_ap = x1(5)/C_ap; % pulmonary arteries
P_C_vp = x1(6)/C_vp; % pulmonary veins

%creation of diodes

if P_lv < P_C_vp
    alpha_lv = 1;
else
    alpha_lv = 0;
end

if P_rv < P_C_vs
    alpha_rv = 1;
else
    alpha_rv = 0;
end

if P_lv > P_C_as
    beta_lv = 1;
else
    beta_lv = 0;
end

if P_rv > P_C_ap
    beta_rv = 1;
else
    beta_rv = 0;
end

%assigning state variables

dx1 = zeros(6,1);

dx1(1) = ((P_C_vp - P_lv)/R_vp)*alpha_lv - ((P_lv - P_C_as)/R_cs)*beta_lv; %V_lv
dx1(2) = ((P_C_vs - P_rv)/R_vp)*alpha_rv - ((P_rv - P_C_ap)/R_cp)*beta_rv; %V_rv
dx1(3) = ((P_lv - P_C_as)/R_cs)*beta_lv - ((P_C_as - P_C_vs)/R_as); %V_C_as
dx1(4) = ((P_C_as - P_C_vs)/R_as) - ((P_C_vs - P_rv)/R_vs)*alpha_rv; %V_C_vs
dx1(5) = ((P_rv - P_C_ap)/R_cp)*beta_rv - ((P_C_ap - P_C_vp)/R_ap); %V_C_ap
dx1(6) = ((P_C_ap - P_C_vp)/R_ap) - ((P_C_vp - P_lv)/R_vp)*alpha_lv; %V_C_vp

end

[Torsten]

I would like to shift this function 16ms forward:

Maybe

tm1 = mod(t1+0.016,tc1);
%impulse fxn
if tm1 < 1.5*T_es
    a1 = 0.5*(1 + sin((pi*tm1/T_es) - (pi/2)));
else
    a1 = 0.5*exp(-(tm1 - 1.5*T_es)/tau);
end

?

[Hannah Goldblatt]

This seems to be working, thank you!