%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% Modelos Standard PK/PD Atracurium %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
% Modelos Standard PK/PD para o Atracurium                                %
%                                                                         %
% Miguel Teixeira, Novembro, 2012                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Necessárias as bases de dados com os parâmetros dos modelos Standard para o Atracurium
% (Contém ficheiro Mat_Ind_Par com Parâmetros Vl, Cl, T1, T2, Ke0, C50, Gamma e Tau)
load ([cd ''])

%% Cálculo dos parâmetros necessários para construir os modelos compartimentais
% Parâmetros necessários para calcular os parâmetros da parte linear.
Lambda1 = -log(0.5)/mat_ind_par(i,3);
Lambda2 = -log(0.5)/mat_ind_par(i,4);

% Parametros da Parte Linear
ke0     = mat_ind_par(i,5);                                                % [min^-1]
k20     = 0;                                                               % [min^-1]
k10     = mat_ind_par(i,2)./mat_ind_par(i,1);                              % [min^-1]
k21     = (Lambda1.*Lambda2)./k10;                                         % [min^-1]
k1e     = ke0*(1./mat_ind_par(i,1));                                       % [min^-1]
k12     = -(-Lambda1*k10-Lambda2*k10+Lambda1*Lambda2+(k10^2))/(-k20+k10);  % [min^-1]

% Parametros da Parte Não-Linear
C50     = mat_ind_par(i,6);                                                % [ug/kg]
Gamma   = mat_ind_par(i,7);                                                % [adimenssional]

%% Construção das matrizes do modelo compartimental
Matriz_A = [-(k12+k10+k1e) k21 0;k12 -(k21+k20) 0;k1e 0 -ke0];
Matriz_B = [1 0 0]';
Matriz_C = [0 0 1];
Matriz_D = 0;

% Discretização das Matrizes
[Matriz_Aj, Matriz_Bj] = c2d(Matriz_Ai, Matriz_Bi, Ts);                    % As matrizes C e D discretizadas ficam iguais às contínuas, logo não é necessária a sua discretização.
                                                                           % Ts é o tempo de amostragem em minutos (Ex: 20s é 1/3 de minuto).

%% Sistema de Espaço de Estados
Estado = Matriz_A*Estado+Matriz_B*U;                                       % Estado são as concentrações nos diferentes compartimentos [ug/kg]. U é o Input de Atracurium [ug/kg/min].
Ce = Matriz_C*Estado+Matriz_D;                                             % Obtenção da Concentração de Efeito Ce [ug/kg]
        
%% Equação de Hill (Parte Não Linear)
NMB = 100/(1+((Ce/C50)^Gamma));                                            % Obtenção do Nível de NMB [%]