Merge branch 'dev' into mit_expression

etfl/analysis/dynamic.py

@@ -361,12 +361,14 @@ def update_medium(t, Xi, Si, dmodel, medium_fun, timestep):
 
     return X,S
 
-def compute_center(dmodel, provided_solution=None):
+def compute_center(dmodel, objective,
+                   provided_solution=None, revert_changes=True):
     """
     Fixes growth to be above computed lower bound, finds chebyshev center,
     resets the model, returns solution data
 
     :param dmodel:
+    :param objective: the radius to maximize
     :return:
     """
 
@@ -386,12 +388,13 @@ def compute_center(dmodel, provided_solution=None):
 
     fix_growth(dmodel, the_solution)
 
+    dmodel.objective = objective #dmodel.chebyshev_radius.radius.variable
     try:
-        dmodel.objective = dmodel.chebyshev_radius.radius.variable
         dmodel.optimize()
         chebyshev_sol = dmodel.solution
 
-        release_growth(dmodel)
+        if revert_changes:
+            release_growth(dmodel)
     except AttributeError: #does not solve
         dmodel.logger.warning('Chebyshev solution at fixed growth solution '
                               'not found - relaxing integers and solving '
@@ -403,8 +406,9 @@ def compute_center(dmodel, provided_solution=None):
         dmodel.optimize()
         chebyshev_sol = dmodel.solution
 
-    dmodel.growth_reaction.lower_bound = prev_lb
-    dmodel.objective = prev_obj
+    if revert_changes:
+        dmodel.growth_reaction.lower_bound = prev_lb
+        dmodel.objective = prev_obj
     return chebyshev_sol
 
 
@@ -474,7 +478,9 @@ def run_dynamic_etfl(model, timestep, tfinal, uptake_fun, medium_fun,
 
     # We want to compute the center under the constraint of the growth at the
     # initial solution.
-    chebyshev_sol = compute_center(dmodel,provided_solution=initial_solution)
+    chebyshev_sol = compute_center(dmodel,
+                                   objective = dmodel.chebyshev_radius.radius.variable,
+                                   provided_solution=initial_solution)
 
     # Only now do we add the dynamic variable constraints.
     add_dynamic_variables_constraints(dmodel, timestep, dynamic_constraints)
@@ -561,4 +567,4 @@ def run_dynamic_etfl(model, timestep, tfinal, uptake_fun, medium_fun,
     return wrap_time_sol(var_solutions, obs_values)
 
 def wrap_time_sol(var_solutions, obs_values):
-    return pd.concat([var_solutions,obs_values], axis=0)
+    return pd.concat([var_solutions,obs_values], axis=0)

etfl/core/allocation.py

@@ -399,12 +399,7 @@ def add_dna_mass_requirement(model, mu_values, dna_rel, gc_ratio,
                                  scaled=True)
     model.add_reactions([dna_formation])
 
-    # In this formulation, we make 1 unit of the whole chromosome with NTPs
-    g = gc_ratio
-    mets = {v: -1 * chromosome_len * (g if k.lower() in 'gc' else 1 - g)
-            for k, v in model.dna_nucleotides.items()}
-    # Don't forget to release ppi (2 ppi per bp)
-    mets[ppi] = 2 * chromosome_len
+    mets = get_dna_synthesis_mets(model, chromosome_len, gc_ratio, ppi)
 
     dna_formation.add_metabolites(mets)
 
@@ -443,6 +438,16 @@ def add_dna_mass_requirement(model, mu_values, dna_rel, gc_ratio,
     model.regenerate_constraints()
 
 
+def get_dna_synthesis_mets(model, chromosome_len, gc_ratio, ppi):
+    # In this formulation, we make 1 unit of the whole chromosome with NTPs
+    g = gc_ratio
+    mets = {v: -1 * chromosome_len * (g if k.lower() in 'gc' else 1 - g)
+            for k, v in model.dna_nucleotides.items()}
+    # Don't forget to release ppi (2 ppi per bp)
+    mets[ppi] = 2 * chromosome_len
+    return mets
+
+
 def define_dna_weight_constraint(model, dna, dna_ggdw, gc_content, chromosome_len):
     # DNA mass (BioPython has g.mol^-1, while we are in mmol)
     ma = molecular_weight('A', seq_type='DNA') / 1000  # g.mol^-1 -> kg.mol^-1 (SI) = g.mmol^-1

etfl/core/equilibrium.py

@@ -0,0 +1,121 @@
+from .genes import ExpressedGene
+from ..optim.variables import RNAPUsage, FreeEnzyme, BinaryActivator, \
+    LinearizationVariable
+from ..optim.constraints import GeneConstraint, EnzymeRatio,\
+    LinearizationConstraint
+from pytfa.optim.reformulation import petersen_linearization
+
+
+class RNAPBindingEquilibrium(GeneConstraint):
+    prefix = 'PBEq_'
+
+
+def linearize_product(model, b, x, queue=False):
+    """
+
+    :param model:
+    :param b: the binary variable
+    :param x: the continuous variable
+    :param queue: whether to queue the variables and constraints made
+    :return:
+    """
+
+    # Linearization step for ga_i * [E]
+    z_name = '__MUL__'.join([b.name, x.name])
+    # Add the variables
+    model_z_u = model.add_variable(kind=LinearizationVariable,
+                                  hook=model,
+                                  id_=z_name,
+                                  lb=0,
+                                  ub=x.ub,
+                                  queue=False)
+
+    big_m = x.ub*10
+
+    z_u, new_constraints = petersen_linearization(b=b, x=x, M=big_m,
+                                                  z=model_z_u)
+
+    # Add the constraints:
+    for cons in new_constraints:
+        model.add_constraint(kind=LinearizationConstraint,
+                             hook=model,
+                             id_=cons.name,
+                             expr=cons.expression,
+                             # expr=new_expression,
+                             ub=cons.ub,
+                             lb=cons.lb,
+                             queue=True)
+
+    model._push_queue()
+    return model_z_u
+
+def add_rnap_binding_equilibrium_constraints(model, the_rnap, Kb):
+    """
+
+    :param model:
+    :param the_rnap:
+    :param Kb:
+    :return:
+    """
+
+    # #1. Find if there is a ratio constraint on the RNAP - this is incompatible
+    # all_ratio_cons = model.get_constraints_of_type(EnzymeRatio)
+    #
+    # for rnap_id in model.rnap:
+    #     try:
+    #         model.remove_constraint(all_ratio_cons.get_by_id(rnap_id))
+    #     except KeyError:
+    #         pass
+
+    #2. Create the linearized bilinear product delta_u*RNAP_F
+
+    RNAP_F = model.get_variables_of_type(FreeEnzyme).get_by_id(the_rnap.id).variable
+    activation_vars = model.get_variables_of_type(BinaryActivator)
+
+    zs = {delta_u.name: linearize_product(model, delta_u, RNAP_F, queue=True)
+          for delta_u in activation_vars}
+
+    # TODO: cleanup
+    dna_cons = model.interpolation_constraint.DNA_interpolation.constraint
+    dna_hat = {k.name:v for k,v in
+               dna_cons.get_linear_coefficients(
+                   activation_vars.list_attr('variable')).items()}
+
+    free_rnap_times_gene_conc = sum(zs[delta_u]*dna_hat[delta_u] for delta_u in zs)
+
+    sorted_dna_conc = sorted(dna_hat.values())
+    dna_hat_resolution = max([x-y for x,y in zip(sorted_dna_conc[1:], sorted_dna_conc[:-1])])
+
+
+    all_rnap_usage = model.get_variables_of_type(RNAPUsage)
+
+    #3. For each gene:
+
+    for the_gene in model.genes:
+        if not isinstance(the_gene, ExpressedGene):
+            continue
+
+        copy_number = the_gene.copy_number
+
+        #3.1. get rnap usage
+
+        RNAP_l = all_rnap_usage.get_by_id(the_gene.id)
+
+        #3.2. create constraint
+
+        expr = Kb * RNAP_l - copy_number * free_rnap_times_gene_conc
+
+        # scaling
+        # expr *= model.dna.scaling_factor
+
+        #3.3 Add it to the model
+
+        model.add_constraint(kind=RNAPBindingEquilibrium,
+                             hook=the_gene,
+                             expr=expr,
+                             ub= copy_number*dna_hat_resolution/2,
+                             lb=-copy_number*dna_hat_resolution/2,
+                             queue=True)
+    model._push_queue()
+    model.regenerate_constraints()
+    model.regenerate_variables()

etfl/core/genes.py

@@ -48,7 +48,7 @@ def __init__(self, id, name, sequence, copy_number=1,
         self._rna = ''
         self._peptide = ''
 
-        self._copy_number = copy_number
+        self._copy_number = int(copy_number)
         self._transcribed_by = transcribed_by
         self._translated_by = translated_by
 
@@ -62,10 +62,10 @@ def copy_number(self, value):
         if value != self._copy_number:
             if self.model is None:
                 # Easy
-                self._copy_number = value
+                self._copy_number = int(value)
             else:
                 # We need to make the model change the RNAP allocation
-                self._copy_number = value
+                self._copy_number = int(value)
                 self.model.edit_gene_copy_number(self.id)
         else:
             # Nothing to do here :)

etfl/core/macromolecule.py

@@ -16,11 +16,12 @@
 
 
 class Macromolecule(Species, ABC):
-    def __init__(self, id=None, kdeg=0, *args, **kwargs):
+    def __init__(self, id=None, kdeg=0, scaling_factor = None, *args, **kwargs):
         Species.__init__(self, id = id, *args, **kwargs)
 
         self.kdeg = kdeg
         self.degradation = None
+        self._scaling_factor = scaling_factor
 
     @abstractmethod
     def init_variable(self, queue=False):
@@ -95,4 +96,7 @@ def molecular_weight(self):
 
     @property
     def scaling_factor(self):
-        return 1/self.molecular_weight
+        # Calculate it only once
+        if self._scaling_factor is None:
+            self._scaling_factor = 1/self.molecular_weight
+        return self._scaling_factor

etfl/core/memodel.py

@@ -331,10 +331,11 @@ def add_dummies(self, nt_ratios, mrna_kdeg, mrna_length, aa_ratios,
         add_interpolation_variables(self)
 
         # Create a dummy gene and override the sequences with input data
+        dummy_sequence = 'N'*mrna_length
         dummy_gene = ExpressedGene(id='dummy_gene',
                                    name='Dummy Gene',
-                                   sequence='')
-        dummy_gene._rna = 'N'*mrna_length
+                                   sequence=dummy_sequence)
+        dummy_gene._rna = dummy_sequence
         dummy_gene._peptide = 'X'*peptide_length
 
         self.add_genes([dummy_gene])
@@ -463,7 +464,7 @@ def express_genes(self, gene_list):
             if isinstance(gene, str):
                 the_gene = self.genes.get_by_id(gene)
             else:
-                the_gene = gene
+                the_gene = self.genes.get_by_id(gene.id)
 
             if not isinstance(the_gene, ExpressedGene):
                 continue
@@ -1338,10 +1339,10 @@ def edit_gene_copy_number(self, gene_id):
 
         rnap_alloc = self.get_constraints_of_type(RNAPAllocation)
         rnap_usage = self.get_variables_of_type(RNAPUsage)
-        if self.id in rnap_alloc and self.id in rnap_usage:
+        if gene_id in rnap_alloc and gene_id in rnap_usage:
             # We need to edit the variable
             # Modify the polymerase constraint
-            cons = rnap_alloc.get_by_id(self.id)
+            cons = rnap_alloc.get_by_id(gene_id)
             new_expr = self._get_rnap_allocation_expression(the_gene)
             cons.change_expr(new_expr)
             self.logger.debug('Changed RNAP allocation for gene {}'.format(gene_id))
@@ -1649,7 +1650,6 @@ def _populate_ribosomes(self):
 
         # 3 -> Add ribosomal capacity constraint
         self.regenerate_variables()
-
 
         for rib_id in self.ribosome.keys():
             # CATCH : This is summing ~1500+ variable objects, and for a reason