function [IP3Amount, IP3PkLoc, IP3Amp, IP3TrueDur]=IP3Dyn_model_TH(x,t)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright: Marsa Taheri and Gregory Handy, 2016
% This code was used to simulate the mathematical model of Astrocyte
% IP3dependent Ca responses in 2 papers submitted in Nov 2016.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Given the IP3 trace and corresponding time, it finds several IP3
%characteristics: total amount (area under curve), peak time, peak
%amplitude, and total IP3 duration (from start to where IP3 reached <0.005).
%Multiple things will need to change in this function if the baseline [IP3]
%is not set at 0.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[IP3Amp, IP3PkLocTemp] = findpeaks(x);
IP3PkLoc = t(IP3PkLocTemp);
IP3AmountTemp = cumsum(x);
IP3Amount = IP3AmountTemp(end)/length(x)*t(end); %This is fine for IP3 amount
%because the baseline is 0. If that changes, then may need to change this (like Ca).
%%Or:
%IP3Amount = trapz(x)/length(x)*t(end);
if isempty(IP3PkLocTemp)==0
t_end = t(IP3PkLocTemp + find(x(IP3PkLocTemp:end)<0.005, 1,'first')  1);
else
t_end =[];
end
t_start = t(find(x>0, 1,'first')); %1st point where IP3>0
if isempty(t_end)==0 %if there is an end to the IP3 dynamics
IP3TrueDur = t_end  t_start;
else
IP3TrueDur = t(end); %assume the IP3 is for the duration of the entire
%simulation (including the time before any IP3 was applied!)
end
end
