/*
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
*/
/* A Scalable Algorithm to Explore the Gibbs Energy Landscape of
Genoma-Scale Metabolic Networks
D. De Martino, M. Figliuzzi, A. De Martino and E. Marinari,
PLoS Computational Biology,
code released May 2012.
*/
// This code corrects thermodynamic unfeasibilities among fluxes
// and chemical potentials in constraint-based models of metabolic networks
#include "iostream"
#include "math.h"
#include "stdlib.h"
#include "time.h"
#include "fstream"
#include "vector"
#include "string"
#include "sstream"
using namespace std;
int sign(double v){
int s = ( v > 0 ) ? 1:-1;
return s;
}
void usage(){
cout << " Usage: -m -p -f " << endl;
}
int main (int argc,char** argv){
int M = 0;
int N = 0;
int ok = 0;
// reading the input files
fstream file;
char *matrice = NULL;
char *potenziali = NULL;
char *flussi = NULL;
int c;
while((c = getopt(argc, argv, "m:p:f:"))!=-1) {
switch(c) {
case 'm':
matrice = optarg;
break;
case 'p':
potenziali = optarg;
break;
case 'f':
flussi = optarg;
break;
default:
cout << "error parsing command line" << endl;
usage();
exit(1);
}
}
file.open(matrice, ios::in);
if(file.is_open()==0){
cout << "the file " << matrice << " doesn't exist" << endl;
exit(1);
}
do{
string papa;
while(file.eof()!=1){
getline(file,papa); // computing the size of the stochiometric matrix
N++;
}
N--;
file.close();
file.open(matrice, ios::in);
double a;
while(file >> a) M++;
M/=N;
file.close();
}while(ok==1);
std::vector react[N];
std::vector metab[M];
std::vector react2[N]; // definition of the stochiometrix matrix, flux and chemical potential vector
std::vector metab2[M];
double flux[N];
double chempot0[M], chempot[M];
do{
file.open(matrice, ios::in); // reading the stochiometric matrix
int i1 = 0;
int j1 = 0;
int y=0;
double pepe;
while(file >> pepe){
if( pepe != 0 ){
react[i1].push_back(j1);
react2[i1].push_back(pepe);
metab[j1].push_back(i1);
metab2[j1].push_back(pepe);
}
j1++;
if( j1 == M ){
j1 = 0;
i1++;
}
}
file.close();
file.open(flussi,ios::in); // reading the flux vector
if(file.is_open()==0){
cout << "the file " << flussi << " doesn't exist" << endl;
exit(1);
}
for(int i = 0; i <= N-1; i++) file >> flux[i];
file.close();
file.open(potenziali,ios::in); // reading the chemical potential vector
if(file.is_open()==0){
cout << "the file " << potenziali << " doesn't exist" << endl;
exit(1);
}
for(int i = 0; i <= M-1; i++) file >> chempot0[i];
file.close();
for(int i = 0; i <= M-1; i++) chempot[i] = chempot0[i];
} while(ok==1);
int count = 0;
int ourMax = 30 * N;
// ourMax is the time when the algorithm starts searching for loops
int unsat[ourMax];
do{
if( count == ourMax ){
// for exceedingly large times search for loops among the unsat reactions
cout << "Exceedingly large times: searching for loops" << endl;
count = 0;
std::vector distinti;
distinti.push_back(unsat[0]);
for(int i = 1; i <= ourMax-1; i++){
int oko = 0;
for(int j = 0; j <= i-1; j++) if(unsat[i] == unsat[j])oko++;
if(oko == 0) distinti.push_back(unsat[i]);
}
std::vector matched;
for(int i = 0; i <= distinti.size()-1; i++){
int r = distinti[i];
// search for distinct and connected reactions among the unsat
int match = 1;
for(int j = 0; j <= react[r].size()-1; j++){
int k = 0;
int okei = 0;
do{
if(k != i){
int r1 = distinti[k];
for(int j1 = 0; j1 <= react[r1].size()-1; j1++){
if(react[r][j] == react[r1][j1]) okei = 1;
}
}
k++;
} while( (k0) ){
cout << "searching among " << matched.size() << " reactions" << endl;
std::vector cycle;
int check = 0;
int lenght = 2; // exaustive search (max lenght 26)
int index[matched.size()];
for(int i = 0; i <= matched.size()-1; i++) index[i] = i;
cout << "checking for loops large " << lenght << endl;
do{
int nodes[lenght];
for(int i = 0; i <= lenght-1; i++) nodes[i] = matched[index[i]];
int metab = 0;
int bell = 0;
do{
for(int i = 0; i <= lenght-1; i++){
int r = nodes[i];
int s = 1;
if(flux[r]<0) s = -1;
for(int j = 0; j <= react[r].size()-1; j++){
if(react[r][j] == metab) bell += s * react2[r][j];
}
}
metab++;
} while( (metabmatched.size()-i){
index[lenght-i-1]++;
for(int k = 0; k <= i-1; k++){
index[lenght-i+k] = index[lenght-i-1+k]+1;
}
}
}
if(index[0]>matched.size()-lenght){
lenght++;
if( (lenght <= matched.size()) && (check == 0) ){
cout << "checking for loops large " << lenght << endl;
}
for(int i = 0; i <= matched.size()-1; i++) index[i] = i;
}
} while( (check == 0) && (lenght= 0) ok = 1;
else{
cout << "reaction "<< min << " has flux v_" << min << "=" << flux[min] << " but the free energy is G_" << min << "=" << -xmin * sign(flux[min]) << endl;
double alpha = 0.1; // alpha is the lenght of the update step
for(int j = 0; j <= react[min].size()-1; j++){
int cip = react[min][j];
cout << "from g_" << cip+1 << "=" << chempot[cip] ;
if(flux[min]>0) chempot[cip] -= alpha * react2[min][j]; // update
else chempot[cip] += alpha * react2[min][j];
cout << " to g_" << cip+1 << "=" << chempot[cip] << endl;
}
}
} while(ok == 0);
cout << "Chemical potentials: "<< endl;
for(int j = 0; j <= M-1; j++) cout << chempot[j] << " ";
// writing a thermodynamically feasible chemical potential and flux vector
cout << endl;
cout << "Fluxes: "<< endl;
for(int j = 0; j <= N-1; j++) cout << flux[j] << " ";
cout << endl;
return 0 ;
}