classdef WallMechanics < handle
% The 'WallMechanics' code contains the following sections of the model:
% The Contraction and Mechanical Models.
% Please refer to the relevient sections in the documentation for
% full information on the equations and variable names.
properties
params
u0
index
n_out
idx_out
enabled
end
methods
function self = WallMechanics(varargin)
self.params = parse_inputs(varargin{:});
self.index = indices();
self.u0 = initial_conditions(self.index);
self.enabled = true(size(self.u0));
[self.idx_out, self.n_out] = output_indices();
end
function [du, varargout] = rhs(self, t, u, Ca_i, R_cGMP2)
% Initalise inputs and parameters
idx = self.index;
p = self.params;
[R, h] = self.shared(t, u);
Mp = u(idx.Mp, :);
AMp = u(idx.AMp, :);
AM = u(idx.AM, :);
%% Contraction Equations
K_1 = p.gamma_cross * Ca_i.^p.n_cross;
K_6 = K_1;
K_2 = 58.1395 * p.k_mlcp_b + 58.1395 * p.k_mlcp_c * R_cGMP2;
K_5 = K_2;
M = 1 - AM - AMp - Mp;
du(idx.Mp, :) = p.wallMech * ( p.K_4 * AMp + K_1 .* M - (K_2 + p.K_3) .* Mp );
du(idx.AMp, :) = p.wallMech * ( p.K_3 * Mp + K_6 .* AM - (p.K_4 + K_5) .* AMp );
du(idx.AM, :) = p.wallMech * ( K_5 .* AMp - (p.K_7 + K_6) .* AM );
% Mechanical Equations
F_r = AMp + AM;
E = p.E_passive + F_r * (p.E_active - p.E_passive);
R_0 = p.R_0_passive + F_r * (p.alpha - 1) * p.R_0_passive;
du(idx.R, :) = p.R_0_passive / p.eta * ( R * p.trans_p ./ h - E .* (R - R_0) ./ R_0);
du = bsxfun(@times, self.enabled, du);
if nargout == 2
Uout = zeros(self.n_out, size(u, 2));
Uout(self.idx_out.M, :) = M;
Uout(self.idx_out.F_r, :) = F_r;
Uout(self.idx_out.R, :) = R;
Uout(self.idx_out.K_2, :) = K_2;
varargout{1} = Uout;
end
end
function [R, h] = shared(self, ~,u)
R = u(self.index.R, :);
h = 0.1 * R;
end
function names = varnames(self)
names = [fieldnames(self.index); fieldnames(self.idx_out)];
end
end
end
function idx = indices()
% Index of parameters needing inital conditions
idx.Mp = 1;
idx.AMp = 2;
idx.AM = 3;
idx.R = 4;
end
function [idx, n] = output_indices()
% Index of all other output parameters
idx.M = 1;
idx.F_r = 2;
idx.R = 3;
idx.K_2 = 4;
n = numel(fieldnames(idx));
end
function params = parse_inputs(varargin)
parser = inputParser();
% Parameter for changing the wall mechanics rate constants
parser.addParameter('wallMech', 1);
% Contraction Equation Constants
parser.addParameter('K_3', 0.4); % s^-1
parser.addParameter('K_4', 0.1); % s^-1
parser.addParameter('K_7', 0.1); % s^-1
parser.addParameter('gamma_cross', 17); %uM^-3 s^-1
parser.addParameter('n_cross', 3); % fraction constant of the phosphorylation crossbridge
% Mechanical Equation Constants
parser.addParameter('eta', 1e4); %Pa s
parser.addParameter('R_0_passive', 20e-6); % m
parser.addParameter('trans_p', 4000); % Pa transmural pressure
parser.addParameter('E_passive', 66e3); % Pa
parser.addParameter('E_active', 233e3); % Pa
parser.addParameter('alpha', 0.6); % [-]
parser.addParameter('k_mlcp_b', 0.0086); % [s^-1]
parser.addParameter('k_mlcp_c', 0.0327); % [s^-1]
parser.parse(varargin{:});
params = parser.Results;
end
function u0 = initial_conditions(idx)
u0 = zeros(length(fieldnames(idx)), 1);
% Inital estimations of parameters from experimental data
u0(idx.Mp) = 0.0842; % [-]
u0(idx.AMp) = 0.0622; % [-]
u0(idx.AM) = 0.2746; % [-]
u0(idx.R) = 22.97e-6; % [um]
end