-1

I have a set of differential equations in C(created from a tool that takes an xml input), in the following form

#ifdef ODEs
dx[0] = Function1(p[67], p[64], p[66], p[65], p[23], x_c[0], x_c[3], p[49]);//
dx[1] = Function2(p[62], p[64], p[66], p[65], p[23], x_c[1], x_c[3], p[40]);//
#endif /* ODEs */

I'm trying to solve a matrix differential equation of the form X' = F(X) + B, in MATLAB. From the above sample code what I have is the content for X' = F(X).

I have created the B matrix in MATLAB. I'm using MinGW in MATLAB to integrate the C code with MATLAB.

I'm facing challenge in converting the set of ODE's in C to a matrix-differential equation form(X' = F(X)), which I can use in MATLAB.

I would like to ask for suggestions on how this can be done.

Edit 1: As an alternative, would it b possible to import the C code with differential equations and later for the matrix in MATLAB?

Edit 2:

#ifdef SIZE_DEFINITIONS
#define N_METABS 13
#define N_ODE_METABS 0
#define N_INDEP_METABS 5
#define N_COMPARTMENTS 1
#define N_GLOBAL_PARAMS 0
#define N_KIN_PARAMS 54
#define N_REACTIONS 11

#define N_ARRAY_SIZE_P  63  // number of parameters
#define N_ARRAY_SIZE_X  5   // number of initials
#define N_ARRAY_SIZE_Y  0   // number of assigned elements
#define N_ARRAY_SIZE_XC 5   // number of x concentration
#define N_ARRAY_SIZE_PC 8   // number of p concentration
#define N_ARRAY_SIZE_YC 0   // number of y concentration
#define N_ARRAY_SIZE_DX 5   // number of ODEs 
#define N_ARRAY_SIZE_CT 0   // number of conserved totals

#endif // SIZE_DEFINITIONS

#ifdef TIME
#define T  <set here a user name for the time variable> 
#endif // TIME

#ifdef NAME_ARRAYS
const char* p_names[] = {"Sucvac", "glycolysis", "phos", "UDP", "ADP", "ATP", "Glcex", "Fruex", "cell", "Vmax1", "Km1Fruex", "Ki1Fru", "Vmax2", "Km2Glcex", "Ki2Glc", "Vmax3", "Km3Glc", "Km3ATP", "Km4Fru", "Ki3G6P", "Ki4F6P", "Vmax4", "Km4Fru", "Km4ATP", "Km3Glc", "Ki3G6P", "Ki4F6P", "Vmax5", "Ki5Fru", "Km5Fru", "Km5ATP", "Ki5ADP", "Vmax6f", "Keq6", "Ki6Suc6P", "Km6F6P", "Ki6Pi", "Ki6UDPGlc", "Km6UDPGlc", "Vmax6r", "Km6UDP", "Km6Suc6P", "Ki6F6P", "Vmax7", "Km7Suc6P", "Vmax8f", "Keq8", "Ki8Fru", "Km8Suc", "Ki8UDP", "Km8UDP", "Vmax8r", "Km8UDPGlc", "Km8Fru", "Ki8Suc", "Vmax9", "Ki9Glc", "Km9Suc", "Ki9Fru", "Vmax10", "Km10F6P", "Vmax11", "Km11Suc",  "" };
const char* x_names[] = {"HexP", "Fru", "Suc", "Glc", "Suc6P",  "" };
const char* y_names[] = { "" };
const char* xc_names[] = {"HexP", "Fru", "Suc", "Glc", "Suc6P",  "" };
const char* pc_names[] = {"Sucvac", "glycolysis", "phos", "UDP", "ADP", "ATP", "Glcex", "Fruex",  "" };
const char* yc_names[] = { "" };
const char* dx_names[] = {"ODE HexP", "ODE Fru", "ODE Suc", "ODE Glc", "ODE Suc6P",  "" };
const char* ct_names[] = { "" };
#endif // NAME_ARRAYS

#ifdef INITIAL
x[0] = 1;   //metabolite 'HexP': reactions
x[1] = 1;   //metabolite 'Fru': reactions
x[2] = 1;   //metabolite 'Suc': reactions
x[3] = 1;   //metabolite 'Glc': reactions
x[4] = 1;   //metabolite 'Suc6P': reactions
#endif /* INITIAL */

#ifdef FIXED
p[0] = 0;   //metabolite 'Sucvac': fixed
p[1] = 0;   //metabolite 'glycolysis': fixed
p[2] = 5.1; //metabolite 'phos': fixed
p[3] = 0.2; //metabolite 'UDP': fixed
p[4] = 0.2; //metabolite 'ADP': fixed
p[5] = 1;   //metabolite 'ATP': fixed
p[6] = 5;   //metabolite 'Glcex': fixed
p[7] = 5;   //metabolite 'Fruex': fixed
p[8] = 1;   //compartment 'cell':fixed
p[9] = 0.286;   //reaction 'v1':  kinetic parameter 'Vmax1'
p[10] = 0.2;    //reaction 'v1':  kinetic parameter 'Km1Fruex'
p[11] = 1;  //reaction 'v1':  kinetic parameter 'Ki1Fru'
p[12] = 0.286;  //reaction 'v2':  kinetic parameter 'Vmax2'
p[13] = 0.2;    //reaction 'v2':  kinetic parameter 'Km2Glcex'
p[14] = 1;  //reaction 'v2':  kinetic parameter 'Ki2Glc'
p[15] = 0.197;  //reaction 'v3':  kinetic parameter 'Vmax3'
p[16] = 0.07;   //reaction 'v3':  kinetic parameter 'Km3Glc'
p[17] = 0.25;   //reaction 'v3':  kinetic parameter 'Km3ATP'
p[18] = 10; //reaction 'v3':  kinetic parameter 'Km4Fru'
p[19] = 0.1;    //reaction 'v3':  kinetic parameter 'Ki3G6P'
p[20] = 10; //reaction 'v3':  kinetic parameter 'Ki4F6P'
p[21] = 0.197;  //reaction 'v4':  kinetic parameter 'Vmax4'
p[22] = 10; //reaction 'v4':  kinetic parameter 'Km4Fru'
p[23] = 0.25;   //reaction 'v4':  kinetic parameter 'Km4ATP'
p[24] = 0.07;   //reaction 'v4':  kinetic parameter 'Km3Glc'
p[25] = 0.1;    //reaction 'v4':  kinetic parameter 'Ki3G6P'
p[26] = 10; //reaction 'v4':  kinetic parameter 'Ki4F6P'
p[27] = 0.164;  //reaction 'v5':  kinetic parameter 'Vmax5'
p[28] = 12; //reaction 'v5':  kinetic parameter 'Ki5Fru'
p[29] = 0.1;    //reaction 'v5':  kinetic parameter 'Km5Fru'
p[30] = 0.085;  //reaction 'v5':  kinetic parameter 'Km5ATP'
p[31] = 2;  //reaction 'v5':  kinetic parameter 'Ki5ADP'
p[32] = 0.379;  //reaction 'v6':  kinetic parameter 'Vmax6f'
p[33] = 10; //reaction 'v6':  kinetic parameter 'Keq6'
p[34] = 0.07;   //reaction 'v6':  kinetic parameter 'Ki6Suc6P'
p[35] = 0.6;    //reaction 'v6':  kinetic parameter 'Km6F6P'
p[36] = 3;  //reaction 'v6':  kinetic parameter 'Ki6Pi'
p[37] = 1.4;    //reaction 'v6':  kinetic parameter 'Ki6UDPGlc'
p[38] = 1.8;    //reaction 'v6':  kinetic parameter 'Km6UDPGlc'
p[39] = 0.2;    //reaction 'v6':  kinetic parameter 'Vmax6r'
p[40] = 0.3;    //reaction 'v6':  kinetic parameter 'Km6UDP'
p[41] = 0.1;    //reaction 'v6':  kinetic parameter 'Km6Suc6P'
p[42] = 0.4;    //reaction 'v6':  kinetic parameter 'Ki6F6P'
p[43] = 0.5;    //reaction 'v7':  kinetic parameter 'Vmax7'
p[44] = 0.1;    //reaction 'v7':  kinetic parameter 'Km7Suc6P'
p[45] = 0.677;  //reaction 'v8':  kinetic parameter 'Vmax8f'
p[46] = 5;  //reaction 'v8':  kinetic parameter 'Keq8'
p[47] = 4;  //reaction 'v8':  kinetic parameter 'Ki8Fru'
p[48] = 50; //reaction 'v8':  kinetic parameter 'Km8Suc'
p[49] = 0.3;    //reaction 'v8':  kinetic parameter 'Ki8UDP'
p[50] = 0.3;    //reaction 'v8':  kinetic parameter 'Km8UDP'
p[51] = 0.3;    //reaction 'v8':  kinetic parameter 'Vmax8r'
p[52] = 0.3;    //reaction 'v8':  kinetic parameter 'Km8UDPGlc'
p[53] = 4;  //reaction 'v8':  kinetic parameter 'Km8Fru'
p[54] = 40; //reaction 'v8':  kinetic parameter 'Ki8Suc'
p[55] = 0.372;  //reaction 'v9':  kinetic parameter 'Vmax9'
p[56] = 15; //reaction 'v9':  kinetic parameter 'Ki9Glc'
p[57] = 10; //reaction 'v9':  kinetic parameter 'Km9Suc'
p[58] = 15; //reaction 'v9':  kinetic parameter 'Ki9Fru'
p[59] = 0.1;    //reaction 'v10':  kinetic parameter 'Vmax10'
p[60] = 0.2;    //reaction 'v10':  kinetic parameter 'Km10F6P'
p[61] = 1;  //reaction 'v11':  kinetic parameter 'Vmax11'
p[62] = 100;    //reaction 'v11':  kinetic parameter 'Km11Suc'
#endif /* FIXED */

#ifdef ASSIGNMENT
x_c[0] = x[0]/p[8]; //concentration of metabolite 'HexP': reactions
x_c[1] = x[1]/p[8]; //concentration of metabolite 'Fru': reactions
x_c[2] = x[2]/p[8]; //concentration of metabolite 'Suc': reactions
x_c[3] = x[3]/p[8]; //concentration of metabolite 'Glc': reactions
x_c[4] = x[4]/p[8]; //concentration of metabolite 'Suc6P': reactions
p_c[0] = p[0]/p[8]; //concentration of metabolite 'Sucvac': fixed
p_c[1] = p[1]/p[8]; //concentration of metabolite 'glycolysis': fixed
p_c[2] = p[2]/p[8]; //concentration of metabolite 'phos': fixed
p_c[3] = p[3]/p[8]; //concentration of metabolite 'UDP': fixed
p_c[4] = p[4]/p[8]; //concentration of metabolite 'ADP': fixed
p_c[5] = p[5]/p[8]; //concentration of metabolite 'ATP': fixed
p_c[6] = p[6]/p[8]; //concentration of metabolite 'Glcex': fixed
p_c[7] = p[7]/p[8]; //concentration of metabolite 'Fruex': fixed
#endif /* ASSIGNMENT */

#ifdef FUNCTIONS_HEADERS
double FunctionForV1(double prod_0, double sub_0, double param_0, double param_1, double param_2); 
double FunctionForV2(double prod_0, double sub_0, double param_0, double param_1, double param_2); 
double FunctionForV3(double sub_0, double modif_0, double sub_1, double prod_0, double param_0, double param_1, double param_2, double param_3, double param_4, double param_5); 
double FunctionForV4(double sub_0, double sub_1, double modif_0, double prod_0, double param_0, double param_1, double param_2, double param_3, double param_4, double param_5); 
double FunctionForV5(double prod_0, double sub_0, double sub_1, double param_0, double param_1, double param_2, double param_3, double param_4); 
double FunctionForV6(double sub_0, double param_0, double param_1, double param_2, double param_3, double param_4, double param_5, double param_6, double param_7, double param_8, double prod_0, double prod_1, double param_9, double param_10, double modif_0); 
double FunctionForV7(double param_0, double sub_0, double param_1); 
double FunctionForV8(double sub_0, double sub_1, double param_0, double param_1, double param_2, double param_3, double param_4, double param_5, double param_6, double param_7, double prod_0, double prod_1, double param_8, double param_9); 
double FunctionForV9(double prod_0, double prod_1, double param_0, double param_1, double param_2, double sub_0, double param_3); 
double FunctionForV10(double sub_0, double param_0, double param_1); 
double FunctionForV11(double param_0, double sub_0, double param_1); 
#endif /* FUNCTIONS_HEADERS */

#ifdef FUNCTIONS
double FunctionForV1(double prod_0, double sub_0, double param_0, double param_1, double param_2)   //Function for v1
{return  param_2*sub_0/(param_1*(1.00000000000000000+prod_0/param_0)+sub_0);} 
double FunctionForV2(double prod_0, double sub_0, double param_0, double param_1, double param_2)   //Function for v2
{return  param_2*sub_0/(param_1*(1.00000000000000000+prod_0/param_0)+sub_0);} 
double FunctionForV3(double sub_0, double modif_0, double sub_1, double prod_0, double param_0, double param_1, double param_2, double param_3, double param_4, double param_5)     //Function for v3
{return  param_5*(sub_1/param_3)*(sub_0/param_2)/((1.00000000000000000+sub_0/param_2)*(1.00000000000000000+sub_1/param_3+modif_0/param_4+0.11300000000000000*prod_0/param_0+0.05750000000000000*prod_0/param_1));} 
double FunctionForV4(double sub_0, double sub_1, double modif_0, double prod_0, double param_0, double param_1, double param_2, double param_3, double param_4, double param_5)     //Function for v4
{return  param_5*(sub_1/param_4)*(sub_0/param_3)/((1.00000000000000000+sub_0/param_3)*(1.00000000000000000+modif_0/param_2+sub_1/param_4+0.11300000000000000*prod_0/param_0+0.05750000000000000*prod_0/param_1));} 
double FunctionForV5(double prod_0, double sub_0, double sub_1, double param_0, double param_1, double param_2, double param_3, double param_4)     //Function for v5
{return  param_4/(1.00000000000000000+sub_1/param_1)*(sub_1/param_3)*(sub_0/param_2)/(1.00000000000000000+sub_1/param_3+sub_0/param_2+sub_1*sub_0/(param_3*param_2)+prod_0/param_0);} 
double FunctionForV6(double sub_0, double param_0, double param_1, double param_2, double param_3, double param_4, double param_5, double param_6, double param_7, double param_8, double prod_0, double prod_1, double param_9, double param_10, double modif_0)   //Function for v6
{return  param_9*(0.05750000000000000*sub_0*0.82310000000000005*sub_0-prod_0*prod_1/param_0)/(0.05750000000000000*sub_0*0.82310000000000005*sub_0*(1.00000000000000000+prod_0/param_3)+param_5*(1.00000000000000000+modif_0/param_2)*(0.82310000000000005*sub_0+param_4)+param_8*0.05750000000000000*sub_0+param_9/(param_10*param_0)*(param_7*prod_0*(1.00000000000000000+0.82310000000000005*sub_0/param_4)+prod_1*(param_6*(1.00000000000000000+param_8*0.05750000000000000*sub_0/(param_4*param_5*(1.00000000000000000+modif_0/param_2)))+prod_0*(1.00000000000000000+0.05750000000000000*sub_0/param_1))));} 
double FunctionForV7(double param_0, double sub_0, double param_1)  //Function for v7
{return  param_1*sub_0/(param_0+sub_0);} 
double FunctionForV8(double sub_0, double sub_1, double param_0, double param_1, double param_2, double param_3, double param_4, double param_5, double param_6, double param_7, double prod_0, double prod_1, double param_8, double param_9)  //Function for v8
{return  (-param_8)*(prod_0*prod_1-sub_0*0.82310000000000005*sub_1/param_0)/(prod_0*prod_1*(1.00000000000000000+sub_0/param_1)+param_5*(prod_1+param_3)+param_6*prod_0+param_8/(param_9*param_0)*(param_7*sub_0*(1.00000000000000000+prod_1/param_3)+0.82310000000000005*sub_1*(param_4*(1.00000000000000000+param_6*prod_0/(param_3*param_5))+sub_0*(1.00000000000000000+prod_0/param_2))));} 
double FunctionForV9(double prod_0, double prod_1, double param_0, double param_1, double param_2, double sub_0, double param_3)    //Function for v9
{return  param_3/(1.00000000000000000+prod_1/param_1)*sub_0/(param_2*(1.00000000000000000+prod_0/param_0)+sub_0);} 
double FunctionForV10(double sub_0, double param_0, double param_1)     //Function for v10
{return  param_1*0.05750000000000000*sub_0/(param_0+0.05750000000000000*sub_0);} 
double FunctionForV11(double param_0, double sub_0, double param_1)     //Function for v11
{return  param_1*sub_0/(param_0+sub_0);} 
#endif /* FUNCTIONS */

#ifdef ODEs
dx[0] = FunctionForV3(p_c[5], x_c[1], x_c[3], x_c[0], p[19], p[20], p[17], p[16], p[18], p[15])*p[8]+FunctionForV4(p_c[5], x_c[1], x_c[3], x_c[0], p[25], p[26], p[24], p[23], p[22], p[21])*p[8]+FunctionForV5(p_c[4], p_c[5], x_c[1], p[31], p[28], p[30], p[29], p[27])*p[8]-2*FunctionForV6(x_c[0], p[33], p[42], p[36], p[34], p[37], p[35], p[41], p[40], p[38], x_c[4], p_c[3], p[32], p[39], p_c[2])*p[8]-FunctionForV8(x_c[1], x_c[0], p[46], p[47], p[54], p[49], p[53], p[48], p[50], p[52], x_c[2], p_c[3], p[45], p[51])*p[8]-FunctionForV10(x_c[0], p[60], p[59])*p[8];   // 
dx[1] = FunctionForV1(x_c[1], p_c[7], p[11], p[10], p[9])*p[8]-FunctionForV4(p_c[5], x_c[1], x_c[3], x_c[0], p[25], p[26], p[24], p[23], p[22], p[21])*p[8]-FunctionForV5(p_c[4], p_c[5], x_c[1], p[31], p[28], p[30], p[29], p[27])*p[8]-FunctionForV8(x_c[1], x_c[0], p[46], p[47], p[54], p[49], p[53], p[48], p[50], p[52], x_c[2], p_c[3], p[45], p[51])*p[8]+FunctionForV9(x_c[1], x_c[3], p[58], p[56], p[57], x_c[2], p[55])*p[8];  // 
dx[2] = FunctionForV7(p[44], x_c[4], p[43])*p[8]+FunctionForV8(x_c[1], x_c[0], p[46], p[47], p[54], p[49], p[53], p[48], p[50], p[52], x_c[2], p_c[3], p[45], p[51])*p[8]-FunctionForV9(x_c[1], x_c[3], p[58], p[56], p[57], x_c[2], p[55])*p[8]-FunctionForV11(p[62], x_c[2], p[61])*p[8]; // 
dx[3] = FunctionForV2(x_c[3], p_c[6], p[14], p[13], p[12])*p[8]-FunctionForV3(p_c[5], x_c[1], x_c[3], x_c[0], p[19], p[20], p[17], p[16], p[18], p[15])*p[8]+FunctionForV9(x_c[1], x_c[3], p[58], p[56], p[57], x_c[2], p[55])*p[8];    // 
dx[4] = FunctionForV6(x_c[0], p[33], p[42], p[36], p[34], p[37], p[35], p[41], p[40], p[38], x_c[4], p_c[3], p[32], p[39], p_c[2])*p[8]-FunctionForV7(p[44], x_c[4], p[43])*p[8];   // 
#endif /* ODEs */
Natasha
  • 1,111
  • 5
  • 28
  • 66
  • "I'm facing challenge in converting" is vague. Post the code that shows what you have tried and as able, point out the area of difficulty. – chux - Reinstate Monica Sep 13 '18 at 15:08
  • When the form is X' = AX, then Function1 and Function2 are linear. You must get the matrices A and B from their definition. – Daniel1000 Sep 13 '18 at 15:33
  • @DanielFrisch Sorry, the notation(AX) that I am using is wrong . I understand it can be used only for linear system. Apologies. function 1 and 2 are non-linear. Would it be right to represent it in the form X' = F(X)? – Natasha Sep 13 '18 at 15:41
  • @chux I have created the xml file and formed the c file out of it. I'm a beginner in C.I don't know how to create the matrix form of differential equation in C. In MATLAB, `[ dx[0] dx[1]] = [function1 function2]` gives the matrix-differential equation. I request for a sample of how to write a similar form in C. – Natasha Sep 13 '18 at 16:31
  • @Natasha if the elements of the array `dx` are supposed to be functions (i.e. `function1`, `function2`), then I think you will want an array of function pointers. Check out [this answer](https://stackoverflow.com/a/252750/3958521) – bigwillydos Sep 13 '18 at 17:40
  • @Natasha [... and formed the c file out of it.](https://stackoverflow.com/questions/52316511/creating-a-matrix-differential-equation-in-c?noredirect=1#comment91583424_52316511) --> post that C file to add useful info to the post. An [MCVE] would be great. – chux - Reinstate Monica Sep 13 '18 at 17:58
  • Minimal ?? I guess that 90% of the code can be safely discarded. Or, better, discard all code and explain your strategy, instead. –  Sep 14 '18 at 07:48
  • @YvesDaoust Sorry about that. It is not a minimal code. But, I took the simplest case that has 5 ODE's. In my real case, I have more than 50 ODE's. This is what I intend to do, for instance, when there are two functions in,say,MATLAB we can pass the ode's to the other function.I'm looking at similar approach. But , since the ode's are in .c and the other part of my code is in MATLAB I am not sure how to proceed. – Natasha Sep 14 '18 at 08:33

1 Answers1

1

You can perfectly translate this C-Code into MATLAB. Note that I increased all indices +1 because MATLAB indices start with one and not with zero as in C. 

clear all


%% DEFINITIONS

% SIZE_DEFINITIONS
N_METABS = 13;
N_ODE_METABS = 0;
N_INDEP_METABS = 5;
N_COMPARTMENTS = 1;
N_GLOBAL_PARAMS = 0;
N_KIN_PARAMS = 54;
N_REACTIONS = 11;

N_ARRAY_SIZE_P  = 63;  % number of parameters
N_ARRAY_SIZE_X  = 5;   % number of initials
N_ARRAY_SIZE_Y  = 0;   % number of assigned elements
N_ARRAY_SIZE_XC = 5;   % number of x concentration
N_ARRAY_SIZE_PC = 8;   % number of p concentration
N_ARRAY_SIZE_YC = 0;   % number of y concentration
N_ARRAY_SIZE_DX = 5;   % number of ODEs 
N_ARRAY_SIZE_CT = 0;   % number of conserved totals

% NAME_ARRAYS
p_names = {'Sucvac', 'glycolysis', 'phos', 'UDP', 'ADP', 'ATP', 'Glcex', 'Fruex', 'cell', 'Vmax1', 'Km1Fruex', 'Ki1Fru', 'Vmax2', 'Km2Glcex', 'Ki2Glc', 'Vmax3', 'Km3Glc', 'Km3ATP', 'Km4Fru', 'Ki3G6P', 'Ki4F6P', 'Vmax4', 'Km4Fru', 'Km4ATP', 'Km3Glc', 'Ki3G6P', 'Ki4F6P', 'Vmax5', 'Ki5Fru', 'Km5Fru', 'Km5ATP', 'Ki5ADP', 'Vmax6f', 'Keq6', 'Ki6Suc6P', 'Km6F6P', 'Ki6Pi', 'Ki6UDPGlc', 'Km6UDPGlc', 'Vmax6r', 'Km6UDP', 'Km6Suc6P', 'Ki6F6P', 'Vmax7', 'Km7Suc6P', 'Vmax8f', 'Keq8', 'Ki8Fru', 'Km8Suc', 'Ki8UDP', 'Km8UDP', 'Vmax8r', 'Km8UDPGlc', 'Km8Fru', 'Ki8Suc', 'Vmax9', 'Ki9Glc', 'Km9Suc', 'Ki9Fru', 'Vmax10', 'Km10F6P', 'Vmax11', 'Km11Suc'};
x_names = {'HexP', 'Fru', 'Suc', 'Glc', 'Suc6P'};
y_names = { '' };
xc_names = {'HexP', 'Fru', 'Suc', 'Glc', 'Suc6P'};
pc_names = {'Sucvac', 'glycolysis', 'phos', 'UDP', 'ADP', 'ATP', 'Glcex', 'Fruex'};
yc_names = { '' };
dx_names = {'ODE HexP', 'ODE Fru', 'ODE Suc', 'ODE Glc', 'ODE Suc6P'};
ct_names = { '' };

% INITIAL
x = NaN(5,1);
x(0+1) = 1;   %metabolite 'HexP': reactions
x(1+1) = 1;   %metabolite 'Fru': reactions
x(2+1) = 1;   %metabolite 'Suc': reactions
x(3+1) = 1;   %metabolite 'Glc': reactions
x(4+1) = 1;   %metabolite 'Suc6P': reactions

% FIXED
p = NaN(63,1);
p(0+1) = 0;   %metabolite 'Sucvac': fixed
p(1+1) = 0;   %metabolite 'glycolysis': fixed
p(2+1) = 5.1; %metabolite 'phos': fixed
p(3+1) = 0.2; %metabolite 'UDP': fixed
p(4+1) = 0.2; %metabolite 'ADP': fixed
p(5+1) = 1;   %metabolite 'ATP': fixed
p(6+1) = 5;   %metabolite 'Glcex': fixed
p(7+1) = 5;   %metabolite 'Fruex': fixed
p(8+1) = 1;   %compartment 'cell':fixed
p(9+1) = 0.286;   %reaction 'v1':  kinetic parameter 'Vmax1'
p(10+1) = 0.2;    %reaction 'v1':  kinetic parameter 'Km1Fruex'
p(11+1) = 1;  %reaction 'v1':  kinetic parameter 'Ki1Fru'
p(12+1) = 0.286;  %reaction 'v2':  kinetic parameter 'Vmax2'
p(13+1) = 0.2;    %reaction 'v2':  kinetic parameter 'Km2Glcex'
p(14+1) = 1;  %reaction 'v2':  kinetic parameter 'Ki2Glc'
p(15+1) = 0.197;  %reaction 'v3':  kinetic parameter 'Vmax3'
p(16+1) = 0.07;   %reaction 'v3':  kinetic parameter 'Km3Glc'
p(17+1) = 0.25;   %reaction 'v3':  kinetic parameter 'Km3ATP'
p(18+1) = 10; %reaction 'v3':  kinetic parameter 'Km4Fru'
p(19+1) = 0.1;    %reaction 'v3':  kinetic parameter 'Ki3G6P'
p(20+1) = 10; %reaction 'v3':  kinetic parameter 'Ki4F6P'
p(21+1) = 0.197;  %reaction 'v4':  kinetic parameter 'Vmax4'
p(22+1) = 10; %reaction 'v4':  kinetic parameter 'Km4Fru'
p(23+1) = 0.25;   %reaction 'v4':  kinetic parameter 'Km4ATP'
p(24+1) = 0.07;   %reaction 'v4':  kinetic parameter 'Km3Glc'
p(25+1) = 0.1;    %reaction 'v4':  kinetic parameter 'Ki3G6P'
p(26+1) = 10; %reaction 'v4':  kinetic parameter 'Ki4F6P'
p(27+1) = 0.164;  %reaction 'v5':  kinetic parameter 'Vmax5'
p(28+1) = 12; %reaction 'v5':  kinetic parameter 'Ki5Fru'
p(29+1) = 0.1;    %reaction 'v5':  kinetic parameter 'Km5Fru'
p(30+1) = 0.085;  %reaction 'v5':  kinetic parameter 'Km5ATP'
p(31+1) = 2;  %reaction 'v5':  kinetic parameter 'Ki5ADP'
p(32+1) = 0.379;  %reaction 'v6':  kinetic parameter 'Vmax6f'
p(33+1) = 10; %reaction 'v6':  kinetic parameter 'Keq6'
p(34+1) = 0.07;   %reaction 'v6':  kinetic parameter 'Ki6Suc6P'
p(35+1) = 0.6;    %reaction 'v6':  kinetic parameter 'Km6F6P'
p(36+1) = 3;  %reaction 'v6':  kinetic parameter 'Ki6Pi'
p(37+1) = 1.4;    %reaction 'v6':  kinetic parameter 'Ki6UDPGlc'
p(38+1) = 1.8;    %reaction 'v6':  kinetic parameter 'Km6UDPGlc'
p(39+1) = 0.2;    %reaction 'v6':  kinetic parameter 'Vmax6r'
p(40+1) = 0.3;    %reaction 'v6':  kinetic parameter 'Km6UDP'
p(41+1) = 0.1;    %reaction 'v6':  kinetic parameter 'Km6Suc6P'
p(42+1) = 0.4;    %reaction 'v6':  kinetic parameter 'Ki6F6P'
p(43+1) = 0.5;    %reaction 'v7':  kinetic parameter 'Vmax7'
p(44+1) = 0.1;    %reaction 'v7':  kinetic parameter 'Km7Suc6P'
p(45+1) = 0.677;  %reaction 'v8':  kinetic parameter 'Vmax8f'
p(46+1) = 5;  %reaction 'v8':  kinetic parameter 'Keq8'
p(47+1) = 4;  %reaction 'v8':  kinetic parameter 'Ki8Fru'
p(48+1) = 50; %reaction 'v8':  kinetic parameter 'Km8Suc'
p(49+1) = 0.3;    %reaction 'v8':  kinetic parameter 'Ki8UDP'
p(50+1) = 0.3;    %reaction 'v8':  kinetic parameter 'Km8UDP'
p(51+1) = 0.3;    %reaction 'v8':  kinetic parameter 'Vmax8r'
p(52+1) = 0.3;    %reaction 'v8':  kinetic parameter 'Km8UDPGlc'
p(53+1) = 4;  %reaction 'v8':  kinetic parameter 'Km8Fru'
p(54+1) = 40; %reaction 'v8':  kinetic parameter 'Ki8Suc'
p(55+1) = 0.372;  %reaction 'v9':  kinetic parameter 'Vmax9'
p(56+1) = 15; %reaction 'v9':  kinetic parameter 'Ki9Glc'
p(57+1) = 10; %reaction 'v9':  kinetic parameter 'Km9Suc'
p(58+1) = 15; %reaction 'v9':  kinetic parameter 'Ki9Fru'
p(59+1) = 0.1;    %reaction 'v10':  kinetic parameter 'Vmax10'
p(60+1) = 0.2;    %reaction 'v10':  kinetic parameter 'Km10F6P'
p(61+1) = 1;  %reaction 'v11':  kinetic parameter 'Vmax11'
p(62+1) = 100;    %reaction 'v11':  kinetic parameter 'Km11Suc'


% ASSIGNMENT
x_c = NaN(5,1);
x_c(0+1) = x(0+1)/p(8+1); %concentration of metabolite 'HexP': reactions
x_c(1+1) = x(1+1)/p(8+1); %concentration of metabolite 'Fru': reactions
x_c(2+1) = x(2+1)/p(8+1); %concentration of metabolite 'Suc': reactions
x_c(3+1) = x(3+1)/p(8+1); %concentration of metabolite 'Glc': reactions
x_c(4+1) = x(4+1)/p(8+1); %concentration of metabolite 'Suc6P': reactions

p_c = NaN(8,1);
p_c(0+1) = p(0+1)/p(8+1); %concentration of metabolite 'Sucvac': fixed
p_c(1+1) = p(1+1)/p(8+1); %concentration of metabolite 'glycolysis': fixed
p_c(2+1) = p(2+1)/p(8+1); %concentration of metabolite 'phos': fixed
p_c(3+1) = p(3+1)/p(8+1); %concentration of metabolite 'UDP': fixed
p_c(4+1) = p(4+1)/p(8+1); %concentration of metabolite 'ADP': fixed
p_c(5+1) = p(5+1)/p(8+1); %concentration of metabolite 'ATP': fixed
p_c(6+1) = p(6+1)/p(8+1); %concentration of metabolite 'Glcex': fixed
p_c(7+1) = p(7+1)/p(8+1); %concentration of metabolite 'Fruex': fixed



% FUNCTIONS
FunctionForV1 = @(prod_0, sub_0, param_0, param_1, param_2) ... %Function for v1
  param_2*sub_0/(param_1*(1.00000000000000000+prod_0/param_0)+sub_0);

FunctionForV2 = @(prod_0, sub_0, param_0, param_1, param_2) ... %Function for v2
  param_2*sub_0/(param_1*(1.00000000000000000+prod_0/param_0)+sub_0);

FunctionForV3 = @(sub_0, modif_0, sub_1, prod_0, param_0, param_1, param_2, param_3, param_4, param_5) ... % Function for v3
  param_5*(sub_1/param_3)*(sub_0/param_2)/((1.00000000000000000+sub_0/param_2)*(1.00000000000000000+sub_1/param_3+modif_0/param_4+0.11300000000000000*prod_0/param_0+0.05750000000000000*prod_0/param_1));    %Function for v3

FunctionForV4 = @(sub_0, sub_1, modif_0, prod_0, param_0, param_1, param_2, param_3, param_4, param_5) ... % Function for v4
  param_5*(sub_1/param_4)*(sub_0/param_3)/((1.00000000000000000+sub_0/param_3)*(1.00000000000000000+modif_0/param_2+sub_1/param_4+0.11300000000000000*prod_0/param_0+0.05750000000000000*prod_0/param_1));    %Function for v4

FunctionForV5 = @(prod_0, sub_0, sub_1, param_0, param_1, param_2, param_3, param_4) ...    %Function for v5
  param_4/(1.00000000000000000+sub_1/param_1)*(sub_1/param_3)*(sub_0/param_2)/(1.00000000000000000+sub_1/param_3+sub_0/param_2+sub_1*sub_0/(param_3*param_2)+prod_0/param_0);

FunctionForV6 = @(sub_0, param_0, param_1, param_2, param_3, param_4, param_5, param_6, param_7, param_8, prod_0, prod_1, param_9, param_10, modif_0) ...  %Function for v6
  param_9*(0.05750000000000000*sub_0*0.82310000000000005*sub_0-prod_0*prod_1/param_0)/(0.05750000000000000*sub_0*0.82310000000000005*sub_0*(1.00000000000000000+prod_0/param_3)+param_5*(1.00000000000000000+modif_0/param_2)*(0.82310000000000005*sub_0+param_4)+param_8*0.05750000000000000*sub_0+param_9/(param_10*param_0)*(param_7*prod_0*(1.00000000000000000+0.82310000000000005*sub_0/param_4)+prod_1*(param_6*(1.00000000000000000+param_8*0.05750000000000000*sub_0/(param_4*param_5*(1.00000000000000000+modif_0/param_2)))+prod_0*(1.00000000000000000+0.05750000000000000*sub_0/param_1))));

FunctionForV7 = @(param_0, sub_0, param_1) ... %Function for v7
  param_1*sub_0/(param_0+sub_0);

FunctionForV8 = @(sub_0, sub_1, param_0, param_1, param_2, param_3, param_4, param_5, param_6, param_7, prod_0, prod_1, param_8, param_9) ... %Function for v8
  (-param_8)*(prod_0*prod_1-sub_0*0.82310000000000005*sub_1/param_0)/(prod_0*prod_1*(1.00000000000000000+sub_0/param_1)+param_5*(prod_1+param_3)+param_6*prod_0+param_8/(param_9*param_0)*(param_7*sub_0*(1.00000000000000000+prod_1/param_3)+0.82310000000000005*sub_1*(param_4*(1.00000000000000000+param_6*prod_0/(param_3*param_5))+sub_0*(1.00000000000000000+prod_0/param_2))));

FunctionForV9 = @(prod_0, prod_1, param_0, param_1, param_2, sub_0, param_3) ...   %Function for v9
  param_3/(1.00000000000000000+prod_1/param_1)*sub_0/(param_2*(1.00000000000000000+prod_0/param_0)+sub_0);

FunctionForV10 = @(sub_0, param_0, param_1) ...    %Function for v10
  param_1*0.05750000000000000*sub_0/(param_0+0.05750000000000000*sub_0);

FunctionForV11 = @(param_0, sub_0, param_1) ...    %Function for v11
  param_1*sub_0/(param_0+sub_0);



% ODE Function
odefun = @(t,x_c) [
  FunctionForV3(p_c(5+1), x_c(1+1), x_c(3+1), x_c(0+1), p(19+1), p(20+1), p(17+1), p(16+1), p(18+1), p(15+1))*p(8+1)+FunctionForV4(p_c(5+1), x_c(1+1), x_c(3+1), x_c(0+1), p(25+1), p(26+1), p(24+1), p(23+1), p(22+1), p(21+1))*p(8+1)+FunctionForV5(p_c(4+1), p_c(5+1), x_c(1+1), p(31+1), p(28+1), p(30+1), p(29+1), p(27+1))*p(8+1)-2*FunctionForV6(x_c(0+1), p(33+1), p(42+1), p(36+1), p(34+1), p(37+1), p(35+1), p(41+1), p(40+1), p(38+1), x_c(4+1), p_c(3+1), p(32+1), p(39+1), p_c(2+1))*p(8+1)-FunctionForV8(x_c(1+1), x_c(0+1), p(46+1), p(47+1), p(54+1), p(49+1), p(53+1), p(48+1), p(50+1), p(52+1), x_c(2+1), p_c(3+1), p(45+1), p(51+1))*p(8+1)-FunctionForV10(x_c(0+1), p(60+1), p(59+1))*p(8+1);   %
  FunctionForV1(x_c(1+1), p_c(7+1), p(11+1), p(10+1), p(9+1))*p(8+1)-FunctionForV4(p_c(5+1), x_c(1+1), x_c(3+1), x_c(0+1), p(25+1), p(26+1), p(24+1), p(23+1), p(22+1), p(21+1))*p(8+1)-FunctionForV5(p_c(4+1), p_c(5+1), x_c(1+1), p(31+1), p(28+1), p(30+1), p(29+1), p(27+1))*p(8+1)-FunctionForV8(x_c(1+1), x_c(0+1), p(46+1), p(47+1), p(54+1), p(49+1), p(53+1), p(48+1), p(50+1), p(52+1), x_c(2+1), p_c(3+1), p(45+1), p(51+1))*p(8+1)+FunctionForV9(x_c(1+1), x_c(3+1), p(58+1), p(56+1), p(57+1), x_c(2+1), p(55+1))*p(8+1);  %
  FunctionForV7(p(44+1), x_c(4+1), p(43+1))*p(8+1)+FunctionForV8(x_c(1+1), x_c(0+1), p(46+1), p(47+1), p(54+1), p(49+1), p(53+1), p(48+1), p(50+1), p(52+1), x_c(2+1), p_c(3+1), p(45+1), p(51+1))*p(8+1)-FunctionForV9(x_c(1+1), x_c(3+1), p(58+1), p(56+1), p(57+1), x_c(2+1), p(55+1))*p(8+1)-FunctionForV11(p(62+1), x_c(2+1), p(61+1))*p(8+1); %
  FunctionForV2(x_c(3+1), p_c(6+1), p(14+1), p(13+1), p(12+1))*p(8+1)-FunctionForV3(p_c(5+1), x_c(1+1), x_c(3+1), x_c(0+1), p(19+1), p(20+1), p(17+1), p(16+1), p(18+1), p(15+1))*p(8+1)+FunctionForV9(x_c(1+1), x_c(3+1), p(58+1), p(56+1), p(57+1), x_c(2+1), p(55+1))*p(8+1);    %
  FunctionForV6(x_c(0+1), p(33+1), p(42+1), p(36+1), p(34+1), p(37+1), p(35+1), p(41+1), p(40+1), p(38+1), x_c(4+1), p_c(3+1), p(32+1), p(39+1), p_c(2+1))*p(8+1)-FunctionForV7(p(44+1), x_c(4+1), p(43+1))*p(8+1);   %
];





%% ODE Solution
% https://de.mathworks.com/help/matlab/math/choose-an-ode-solver.html

tspan = [0,50];

[t,x_c_result] = ode45(odefun, tspan, x_c);



%% Plot Results

name = 'ChemicalAnalysis';

fig = figure(2943934);
set(fig, 'Color','white', 'NumberTitle','off', 'Name',name)
clf(fig)
ax = axes(fig);
set(ax, 'XGrid','on', 'YGrid','on', 'XMinorGrid','on', 'YMinorGrid','on', 'NextPlot','add')
xlabel(ax, 'Time in s')
ylabel(ax, 'Concentrations in ?')

ph = plot(ax, t, x_c_result);
set(ph, 'LineWidth',2)

for iLine = 1:length(ph)
  ph(iLine).DisplayName = xc_names{iLine};
end

lg = legend(ax);
lg.Location = 'EastOutside';

In order to save this figure as PDF or PNG, use the export_fig function.

Daniel1000
  • 779
  • 4
  • 11
  • Thanks a lot.Using, NaN(5,1) gives an array of NaN. Can we use zeros instead? From the code that you have posted, the functions are defined in `odefun` using function handle.While looking at examples of matrix differential equations,I came across the `dsolve` option.To use dsolve, I'm using a 5*1 matrix with the functions present in `odefun` with an additional term(`odefun+eye(5)*x_c(1:5)`).`dsolve` requires symbolic representation of all the variables if one can to solve the differential equation using `dsolve`.Is there a simpler way to solve the matrix differential equations?Any suggestions? – Natasha Sep 14 '18 at 15:57
  • You can also use `x=zeros(5,1)`, but it doesn't matter since the values are overridden with ones in the following lines anyway. I did that just to ensure that `p` has the correct size. – Daniel1000 Sep 17 '18 at 07:27
  • With `ode45`, I solve the equation system numerically, and this is what's typically done with these kinds of differential equations. They are highly nonlinear, so I suppose there are analytic solutions only for special simplified cases. And even if there was a solution it would be so complex that it doesn't help you anyway. – Daniel1000 Sep 17 '18 at 07:32
  • If you still want to try with `dsolve`, then you must translate the problem again into its special syntax as it is well documented. Start with some simpler examples that you can also solve by hand and get used to the syntax. But still, I don't think you want or need analytic solutions. Just work with the numerical results from `ode45`. – Daniel1000 Sep 17 '18 at 07:34
  • Thanks a lot. I just realised `dsolve` can be used to obtain the analytical solutions only. In the code that you posted above, `odefun = @(t,x_c) [FunctionForV3;FunctionForV1;FunctionForV7;FunctionForV2;FunctionForV6]` is a vector. Can this function, which contains a vector, be replaced with a matrix? – Natasha Sep 17 '18 at 08:43
  • By matrix , I mean a matrix vector product.Example `odefun=@(t,x_c) [FunctionForV3;FunctionForV1;FunctionForV7;FunctionForV2;FunctionForV6]+eye(5)*x_c(1:5)` .I tried this, but didn't succeed in obtaining the result. – Natasha Sep 17 '18 at 10:20
  • Defining Function = `[FunctionForV3;FunctionForV1;FunctionForV7;FunctionForV2;FunctionForV6]+eye(5)*x_c(1:5)` then using `odefun=@(t,x_c) Function` worked. Thanks a lot – Natasha Sep 17 '18 at 12:46
  • Well, you could have just multiplied the vector entries with the elements of your x_c: `x_c(1)*FunctionForV3(p_c(5+1), x_c(1...` – Daniel1000 Nov 02 '18 at 11:27