implement ODE for idnlgrey parameter estimation

Purpose: estimation of parameters A,B,C in heart rate model by fitting on experimental data. Optimum for current dataset is already known but I'd like to implement it in Matlab.

smax = 185; %max HR smin = 40; %min HR D = 146; % for this dataset

Equation: s'(t) = A*[s(t)-smin]^B * [smax - s(t)]^C * [D - s(t)]^E

where s(t) = heart rate (bpm), smin<s(t)<smax (can be deducted from the equation: the closer s(t) approaches smax or smin, the smaller s'(t) will be)

Normalisation with: sn(t) = (s(t)-smin)/(smax - smin); Dn = (D-smin)/(smax-smin)

Normalised equation: sn'(t) = An*[sn(t)]^B * [1 - sn(t)]^C * [Dn - sn(t)]^E

Initial values: An = 0.3; % realistic: 0<A<1, optimum: 0.45<A<0.65 B = 1.4; % realistic: 0<B<2.5, optimum: 1.55<B<1.75 C = 1.1; % realistic: 0<C<2.5, optimum: 1.60<C<1.80 E = 1;

Dataset is in attachment (data.mat contains 'y' = 1x960 double Heart Rate data, sample frequency = 0.4479 since length of datacollection is around 430 secs.)

I'd like to fit the ODE (normalized) on the data and get A,B and C parameter estimation values. To fit the normalized ODE outcome to the data y, dataset y also needs to be normalized by: y = (y-smin)/(smax-smin).

I tried already following:

data = iddata(y,[],0.4479,'Name','HeartRate');
data.OutputName = 'Heart Rate';
data.OutputUnit = 'bpm';
data.Tstart = 0;
data.TimeUnit = 's';

A = 0.3; % realistic: 0<A<1, optimum: 0.45<A<0.65
B = 1.4; % realistic: 0<B<2.5, optimum: 1.55<B<1.75
C = 1.1; % realistic: 0<C<2.5, optimum: 1.60<C<1.80
E = 1;

smax = 185; %max HR
smin = 40; %min HR
D = mean(y(length(y)-40):y(end));

order         = [1 0 1]; % is this OK???
parameters    = {A,B,C,E,D,smin,smax};
initial_state = [48];
Ts            = 0;
nonlinear_model = idnlgrey('HR_ode',order,parameters,initial_state,Ts);

setpar(nonlinear_model,'Fixed',{false false false true true true true}); % since only A, B and C can vary

nonlinear_model = nlgreyest(data,nonlinear_model,'Display','Full');

I think something is wrong with my 'HR_ode' implementation... :

function [dHR] = HR_ode(t,sin,~,A,B,C,E,Din,smin,smax,varargin)

Anorm = A*(smax-smin)^(B+C+E-1);
snorm = (sin-smin)/(smax-smin);
Dnorm = Din/(smax-smin);

dsnorm = Anorm*((snorm)^B)*((1-snorm)^C)*((Dnorm-snorm)^E);

dHR = dsnorm*(smax-smin)+smin;

end

I spent much time to explain the issue as good as possible. Hopefully someone would be helpful in finding a solution! If you need more information, please let me know.

Thank you!

Regards, Danny