sampling using the maximum entropy probability distribution

[hariszaf]

Up to now, dingo supports sampling using the uniform, the
multivariate exponential and the multivariate Gaussian distributions.

A quite interesting addition from the biologist-point-of-view would be
to enable sampling using the maximum entropy probability distribution.

You may see more about why this would be the case under this study.

"In addition to accounting for fluctuations, the maximum entropy construction provides a principled interpolation between two extremal regimes of metabolic network function. In the “uniform” (no-optimization) limit, no control is exerted over metabolic fluxes: they are selected at random as long as they are permitted by stoichiometry, resulting in broad yet non-trivial flux distributions that support a small, non-zero growth rate. In
the FBA limit, fluxes are controlled precisely to maximize the growth rate, with zero fluctuations. "

For more about max entropy distribution:
https://journals.aps.org/pr/abstract/10.1103/PhysRev.106.620

[hariszaf]

@TolisChal @vissarion we might find useful the following publications:
http://proceedings.mlr.press/v99/straszak19a/straszak19a.pdf
https://journals-aps-org.kuleuven.e-bronnen.be/pre/pdf/10.1103/PhysRevE.96.060401

and probably the most helpful:
https://doi.org/10.1016/j.heliyon.2018.e00596

[hariszaf]

hi again!

Based on the
code trying to figure out what's going on on Statistical mechanics for metabolic networks during steady state growth but also based on their supplementary
I was wondering if what the authors do during the sampling step is simply to adjust the sampling step size based on the flux maximization concept.

The code is not that clean but is it just the // flux maximization part
and you can now generate samples with a distribution that approximates maximum entropy?

void ellipsHR(){                    // hit and run with ellipsoid
    
    for(int i=0;i<=r-1;i++) previous[i]=flux[i];

    for(int yes=0;yes<=r-1;yes++){            // sweep over the axis directions
       
       long double   tp,tm;
       findextrema(tp,tm,yes);


       long double   t=0;
       int count=0;
       do { 
         long double c = cb[yes];
         c*=BETA;
         long double ranvar = casual();
         if(fabs(c)<EPSILON) t = tm+(tp-tm)*ranvar;
         else{                                  //  flux maximization

            if(fabs(c*(tp-tm))<MAXIMUM_EXP) t = tm + log( 1 + ranvar*(exp(c*(tp-tm))-1))/c;

            else{
                    if(c<0) t = tm+log(1-ranvar)/c;
                    else    t = tp+log(ranvar)/c;
             }
        
           }
           count++;  
        } while ( (t-tm<EPSILON || tp-t<EPSILON) && count<10 );

   }

    if(count==10) t=0;
  
          for(int i=0;i<=r-1;i++) flux[i] += t*versori[yes][i];

            int oko=1;
                    for(int i=0;i<=N-1;i++){
                        double flux1=0;
                        if(matrice[i].size()>0) for(int j=0;j<=matrice[i].size()-1;j++) flux1+=matrice2[i][j]*flux[matrice[i][j]];
                                        
                                        if(flux1<boundmin[i] || flux1>boundmax[i])oko=0;
                                            } 

                    if(oko==1)  for(int i=0;i<=r-1;i++) previous[i]=flux[i];
                    else for(int i=0;i<=r-1;i++) flux[i]=previous[i];

}