create_input_files.m

modelOBJ = sbmlimport('cAMP_model_rev19.xml');

nreactions = length(modelOBJ.Reactions);
nspecies = length(modelOBJ.Species);
fid = fopen('cvdx_f.m', 'w');
fidp = fopen('initial_parameters.m', 'w');

fprintf(fid, '%s\n\n\n', 'function yd = cvdx_f(t, y0)');
fprintf(fid, '%s\n\n\n', 'global data spflag timing newISOconc ;');
% print species initial conditions  to derviative and parameter file
for i=1:nspecies
    strnum = int2str(i);
    str = strcat(modelOBJ.Species(i).name, ' =', 'y0(',strnum,')');
    if i==2
    fprintf(fid,'%s\n','if t >= timing & spflag == 1');
    fprintf(fid,'%s;\n',' ISO = newISOconc');
    fprintf(fid,'%s\n','else');
    fprintf(fid,'%s;\n',str);
    fprintf(fid,'%s\n\n','end');
    else
    fprintf(fid,'%s;\n',str);
    end
    
    val = num2str(modelOBJ.Species(i).InitialAmount);
    str = strcat('y0(',strnum,')',' =', val, ' %',modelOBJ.Species(i).Name );
    fprintf(fidp,'%s;\n',str);
end
%     fprintf(fidp,'y0 = transpose(y0);');

fprintf(fidp,'\n');
fprintf(fidp,'\n');
fprintf(fid,'\n');
fprintf(fid,'\n');

%print parameter initial conditions to derivative and parameter file
% get parameters from object file
count = 0;
params = sbioselect(modelOBJ, 'Type', 'parameter');
for i=1:length(params);
        count = count + 1;
        param(i).name = params(i).Name;
        param(i).value = params(i).Value;
        param(i).units = params(i).ValueUnits;
end
count

for i=1:count
    strnum = int2str(i);
    str = strcat(param(i).name, ' =', ' data(', strnum,')');
    fprintf(fid,'%s;\n',str);
    
    val = num2str(param(i).value);
    str = strcat('data(', strnum,')',' =',val,' %',param(i).name);
    fprintf(fidp,'%s;\n',str);
end

fprintf(fid,'\n');
fprintf(fid,'\n');

%get the differential equations from modelOBJ and print to the derivative
%file

for i=1:nspecies
    strnum = int2str(i);
    str = strcat('yd(',strnum,') = ');
    if(modelOBJ.species(i).ConstantAmount == 0)
        for j=1:nreactions
        nreac = length(modelOBJ.Reactions(j).reactants);
        nprod = length(modelOBJ.Reactions(j).products);
            for k=1:nreac
                if(strcmp(modelOBJ.Species(i).name, modelOBJ.Reactions(j).reactants(k).name) == 1)
            str = strcat(str,' -(',modelOBJ.Reactions(j).reactionRate,')');
                end
            end
       
            for k=1:nprod
                if(strcmp(modelOBJ.Species(i).name, modelOBJ.Reactions(j).product(k).name) == 1)
            str = strcat(str,' +(', modelOBJ.Reactions(j).reactionRate,')');
                end
            end
        end 
    else
       str = strcat(str,' 0.0');
    end
    fprintf(fid,'%s;\n',str);
end

fprintf(fid,'\n\n\n');
fprintf(fid,'yd = transpose(yd);\n');

fclose(fid);
fclose(fidp);