- Author:
- WeiweiAi <wai484@aucklanduni.ac.nz>
- Date:
- 2022-08-15 10:06:08+12:00
- Desc:
- Add Ce_BG for chemical species;
Change C_BG to model the electrical capacitor;
Add the first version of build BG script based on stoichiometry matrices.
- Permanent Source URI:
- https://staging.physiomeproject.org/workspace/6bc/rawfile/a1dd1872c24d859834e27b7193a2317e7a4af835/BG_fit/buildBG_matrix.m
%% Read stoichiometric matrices
function buildBG_matrix(fmaxtrixf,fmaxtrixr,modelName,modelPath)
indexN=3;
opts = detectImportOptions(fmaxtrixf);
opts.SelectedVariableNames = [indexN:length(opts.VariableNames)];
N_f = readmatrix(fmaxtrixf,opts);
kappa_name=opts.SelectedVariableNames';
opts.SelectedVariableNames = [indexN-1];
Kname=readmatrix(fmaxtrixf,opts);
opts.SelectedVariableNames = [indexN-2];
Cname=readmatrix(fmaxtrixf,opts); % Component names
opts = detectImportOptions(fmaxtrixr);
opts.SelectedVariableNames = [3:length(opts.VariableNames)];
N_r = readmatrix(fmaxtrixr,opts);
%% Write to CellML
blank='';
% Model
modelfile=[modelPath modelName '.txt'];
fileID = fopen(modelfile,'w');
fprintf(fileID,'def model %s as\n',modelName);
% Unit
unitimport =["J_per_mol","fmol","fmol_per_sec","per_fmol"];
blanks=strcat('',strings(1,length(unitimport)));
unit_all=[blanks;unitimport;unitimport];
fprintf(fileID,'%3s def import using "../cellLib/BG/units_BG.cellml" for\n',blank);
fprintf(fileID,'%7s unit %s using unit %s;\n',unit_all(:));
fprintf(fileID,'%3s enddef;\n',blank);
fprintf(fileID,'\n');
% Import predefined models
Re=strcat('Re',cell(length(kappa_name),1));
importComp=[Cname; Re];
newComp=[Kname;kappa_name];
for i=1:length(importComp)
fprintf(fileID,'%3s def import using "../cellLib/BG/%s_BG.cellml" for\n',blank,importComp{i});
fprintf(fileID,'%7s comp %s using comp %s_BG;\n',blank,newComp{i},importComp{i});
fprintf(fileID,'%3s enddef;\n',blank);
fprintf(fileID,'\n');
end
% Define the component
fprintf(fileID,'%3s def comp %s as\n',blank,modelName);
fprintf(fileID,'%7s var T: kelvin {pub: in, priv: out};\n',blank);
fprintf(fileID,'%7s var t: second {pub: in, priv: out};\n',blank);
% paras
for i=1:length(newComp)
ctype=importComp{i};
switch ctype
case 'C'
fprintf(fileID,'%7s var q_%s_init: fmol {pub: in, priv: out};\n',blank,newComp{i});
fprintf(fileID,'%7s var K_%s: per_fmol {pub: in, priv: out};\n',blank,newComp{i});
fprintf(fileID,'%7s var q_%s: fmol {pub: out, priv: in};\n',blank,newComp{i});
fprintf(fileID,'%7s var mu_%s: J_per_mol {pub: out, priv: in};\n',blank,newComp{i});
fprintf(fileID,'%7s var v_%s: fmol_per_sec {pub: out,priv: out};\n',blank,newComp{i});
case 'Se'
fprintf(fileID,'%7s var q_%s_init: fmol {pub: in, priv: out};\n',blank,newComp{i});
fprintf(fileID,'%7s var K_%s: per_fmol {pub: in, priv: out};\n',blank,newComp{i});
% fprintf(fileID,'%7s var q_%s: fmol {pub: out, priv: in};\n',blank,newComp(i));
fprintf(fileID,'%7s var mu_%s: J_per_mol {pub: out, priv: in};\n',blank,newComp{i});
fprintf(fileID,'%7s var v_%s: fmol_per_sec {pub: out,priv: out};\n',blank,newComp{i});
case 'Re'
fprintf(fileID,'%7s var kappa_%s: fmol_per_sec {pub: in, priv: out};\n',blank,newComp{i});
fprintf(fileID,'%7s var v_%s: fmol_per_sec {pub: out,priv: in};\n',blank,newComp{i});
fprintf(fileID,'%7s var mu_%s_f: J_per_mol {pub: out,priv: out};\n',blank,newComp{i});
fprintf(fileID,'%7s var mu_%s_r: J_per_mol {pub: out,priv: out};\n',blank,newComp{i});
otherwise
disp('The component is not defined')
return
end
end
fprintf(fileID,'\n');
% Equations for conservation laws
for i=1:length(Kname)
findex=find(N_f(i,:));
if ~isempty(findex)
fRe=convertCharsToStrings(kappa_name(findex));
formatSpec = "-%d{dimensionless}*v_%s";
sfRe = strjoin(compose(formatSpec,N_f(i,findex),fRe)');
else
sfRe='';
end
rindex=find(N_r(i,:));
if ~isempty(rindex)
rRe=convertCharsToStrings(kappa_name(rindex));
formatSpec = "+%d{dimensionless}*v_%s";
srRe = strjoin(compose(formatSpec,N_r(i,rindex),rRe)');
else
srRe='';
end
fprintf(fileID,'%7s v_%s = %s%s;\n',blank,Kname{i},sfRe,srRe);
end
for i=1:length(kappa_name)
findex=find(N_f(:,i));
if ~isempty(findex)
fSpecie=convertCharsToStrings(Kname(findex));
formatSpec = "+%d{dimensionless}*mu_%s";
sfSpecie = strjoin(compose(formatSpec,N_f(findex,i),fSpecie)');
fprintf(fileID,'%7s mu_%s_f = %s;\n',blank,kappa_name{i},sfSpecie);
end
rindex=find(N_r(:,i));
if ~isempty(rindex)
rSpecie=convertCharsToStrings(Kname(rindex));
formatSpec = "+%d{dimensionless}*mu_%s";
srSpecie = strjoin(compose(formatSpec,N_r(rindex,i),rSpecie)');
fprintf(fileID,'%7s mu_%s_r = %s;\n',blank,kappa_name{i},srSpecie);
end
end
fprintf(fileID,'%3s enddef;\n',blank);
fprintf(fileID,'\n');
% Group
fprintf(fileID,'%3s def group as encapsulation for\n',blank);
fprintf(fileID,'%7s comp %s incl\n',blank,modelName);
nnewComp=length(newComp);
blanks=strcat('',strings(1,nnewComp));
data_all=[blanks;newComp'];
fprintf(fileID,'%11s comp %s;\n',data_all);
fprintf(fileID,'%7s endcomp;\n',blank);
fprintf(fileID,'%3s enddef;\n',blank);
fprintf(fileID,'\n');
% Mapping
for i=1:length(newComp)
fprintf(fileID,'%3s def map between %s and %s for\n',blank,modelName,newComp{i});
ctype=importComp{i};
switch ctype
case 'C'
varsc=["q","K","v","mu"];
formatSpec = "%s_%s";
names=strcat(newComp{i},strings(length(varsc),1));
varsn=compose(formatSpec,varsc',names);
vars2=["T","t","q_init",varsc];
newq=sprintf("q_%s_init",newComp{i});
vars1=["T","t",newq,varsn'];
case 'Se'
varsc=["K","v","mu"];
formatSpec = "%s_%s";
names=strcat(newComp{i},strings(length(varsc),1));
varsn=compose(formatSpec,varsc',names);
vars2=["T","t","q_init",varsc];
newq=sprintf("q_%s_init",newComp{i});
vars1=["T","t",newq,varsn'];
case 'Re'
vars2=["T","kappa","A_f","A_r","v"];
nkappa=sprintf("kappa_%s",newComp{i});
nmuf=sprintf("mu_%s_f",newComp{i});
nmur=sprintf("mu_%s_r",newComp{i});
nv=sprintf("v_%s",newComp{i});
vars1=["T",nkappa,nmuf,nmur,nv];
otherwise
disp('The component is not defined')
return
end
nvar=length(vars1);
blanks=strcat('',strings(1,nvar));
data_all=[blanks;vars1;vars2];
fprintf(fileID,'%7s vars %s and %s;\n',data_all(:));
fprintf(fileID,'%3s enddef;\n',blank);
fprintf(fileID,'\n');
end
fprintf(fileID,'enddef;\n');
fprintf(fileID,'\n');
fclose(fileID);
end
