Question

У меня есть набор дифференциальных уравнений в C (созданный из инструмента, который принимает ввод XML), в следующем виде

#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 */

Я пытаюсь решить матричное дифференциальное уравнение видаX '= F (X) + B, в MATLAB.Из приведенного выше примера кода у меня есть содержимое для X '= F (X).

Я создал матрицу B в MATLAB.Я использую MinGW в MATLAB для интеграции кода C с MATLAB.

Я сталкиваюсь с проблемой при преобразовании набора ОДУ в C в форму матричного дифференциального уравнения (X '= F (X)), которую я могу использовать в MATLAB.

Iхотел бы спросить предложения о том, как это можно сделать.

Редактировать 1: В качестве альтернативы можно ли импортировать код C с дифференциальными уравнениями и позже для матрицы в MATLAB?

Редактировать 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 */

Daniel1000 · Answer 1 · 14 сентября 2018

Вы можете прекрасно перевести этот C-код в MATLAB. Обратите внимание, что я увеличил все индексы +1, потому что индексы MATLAB начинаются с единицы, а не с нуля, как в 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';

Чтобы сохранить этот рисунок в формате PDF или PNG, используйте функцию export_fig .

Создание матричного дифференциального уравнения в C

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Создание матричного дифференциального уравнения в C

Пожалуйста, войдите или зарегистрируйтесь чтобы ответить на этот вопрос.

1 Ответ

Пожалуйста, войдите или зарегистрируйтесь что бы добавить комментарий.

Похожие темы