function [optFluxes,yield] =  FBA()

    %Stoichiometric matrix
    S = [1,0,-1,-1,0,0,0,-1;0,-1,1,0,0,0,0,-1;0,0,0,1,1,-1,0,0;0,1,0,0,-1,1,1,-1];
    
    %right hand side vector
    b_eq=zeros(4,1);
    
    %lower and upper bounds
    b_lb = zeros(8,1);
    b_ub = 10*ones(8,1);
    
    %objective function
    f = zeros(1,8);
    f(8) = -1;
    
    %zero A->B
    %b_ub(3) = 0;
    
    %zero C -> D
    %b_ub(6) = 0;
    
    %optimization
    [optFluxes,yield] = linprog(f,[],[],S,b_eq,b_lb,b_ub);
    
    