Location: BG_GiProtein @ 9b5f31178365 / MATLAB / WRONG_A_linearised_findNull.m

Author:
Shelley Fong <s.fong@auckland.ac.nz>
Date:
2021-07-22 12:32:44+12:00
Desc:
Initial commit
Permanent Source URI:
https://models.cellml.org/workspace/6f9/rawfile/9b5f31178365f22c5a7e7c2f23cfbc64e674048f/MATLAB/WRONG_A_linearised_findNull.m
Compare... Download Source 

% given dq/dt = A.q + c where A is a square matrix, use null(A) to find ICs
% of q. 
% A is constructed from linearising equations through lm, and
% contains BG parameters 
% kappa is included in AL: no longer kept separate in c

% 20 July WRONG: can't linearise the highly nonlinear system, so can't find
% A, and can't find nullspace. 

% my linearisation here assumed that ln(A+B) = ln(A) + ln(B)


% Kvec = [K.L, K.R, K.G, K.LR, K.a_GTP, K.beta_gamma, K.a_GDP];
current_dir = [cd, '\BG_GPCR_sneyd\MATLAB'];
Idx_backslash = find(current_dir == filesep);
data_dir =  [current_dir filesep 'data' filesep];
output_dir = [current_dir filesep 'output' filesep];

fID = fopen([data_dir 'Kname.txt'], 'r');
Kname = textscan(fID,'%s');
fclose(fID);

params = load([output_dir 'all_params.mat']);
Knamed = params.Knamed;
kappa = params.kappa_named;
N = length(Kname{1});

Kvec = zeros(1,N); % [K.L, K.R, K.G, K.LR, K.a_GTP, K.beta_gamma, K.a_GDP];
for ik=1:N
    Kvec(ik) = Knamed.(Kname{1}{ik});
end

A(1,:) = [-1,-1,0,1,0,0,0];
A(2,:) = [-1,-1,0,1,0,0,0];
A(3,:) = [0,0,-1,-1,1,1,1];
A(4,:) = [1,1,0,-1,1,1,0];
A(5,:) = [0,0,1,-1,-1,-1,1];
A(6,:) = [0,0,1,1,-1,-1,-1];
A(7,:) = [0,0,1,0,1,-1,-2];

% multiply by kappa
A(1,:) = A(1,:) .* kappa.bind;
A(2,:) = A(2,:) .* kappa.bind;
A(3,:) = A(3,:) .* [0,0,kappa.reassoc, kappa.act, kappa.act, kappa.reassoc+kappa.act, kappa.reassoc];
A(4,:) = A(4,:) .* [kappa.bind, kappa.bind, 0, kappa.bind+kappa.act, kappa.act, kappa.act, 0];
A(5,:) = A(5,:) .* [0,0,kappa.act, kappa.act, kappa.act+kappa.hyd, kappa.act, kappa.hyd];
A(6,:) = A(6,:) .* [0,0,kappa.reassoc, kappa.act, kappa.act, kappa.act+kappa.reassoc, kappa.reassoc];
A(7,:) = A(7,:) .* [0,0,kappa.reassoc, 0, kappa.hyd, kappa.reassoc, kappa.hyd+kappa.reassoc];
% multiply by K
for row = 1:length(A)
    A(row,:) = A(row,:).*Kvec;
% %     A(row,:) = exp(A(row,:).*Kvec); % log untransform - most wrong of 3 options.
%     A(row,:) = log(A(row,:).*Kvec); % log transform.
end

Z = null(A);

save([output_dir 'q_initial_conditionsZ.mat'],'Z'); %

if true
    for j=1:N
        fprintf('var q_%s: fmol {init: %g};\n',Kname{1}{j}, Z(j,1));
    end
end


save([output_dir 'A.mat'],'Z','K'); %