Bioscrape and COBRA in a single deterministic cell results in infeasible solution if total_time > 360

Hi folks,

In section 4.1 of Multi-Paradigm-Composites

Setting total_time = 360 works just fine:

from bioscrape_cobra.bioscrape_cobra_deterministic import BioscrapeCOBRAdeterministic
from bioscrape_cobra.simulate import get_bioscrape_cobra_config, simulate_bioscrape_cobra, plot_variables_list_deterministic
from bioscrape_cobra.plot import plot_single, config_single_cell_bioscrape_cobra_topology
from vivarium.plots.topology import plot_topology

# Show full topology here
#TODO wrap this so I can get the Composite before simulating
#bsc_composite_det = BioscrapeCOBRAdeterministic(get_bioscrape_cobra_config())

biocobra_data_det, comp0_det = simulate_bioscrape_cobra(
    initial_glucose=10, #mM
    initial_lactose=200, #mM
    total_time = 360, 
    output_type='timeseries')

But setting total_time=361 results in an infeasible solution, such that multiplying the flux solution None by the flux_scaling parameter results in a TypeError

biocobra_data_det, comp0_det = simulate_bioscrape_cobra(
    initial_glucose=10, #mM
    initial_lactose=200, #mM
    total_time = 361, 
    output_type='timeseries')

Using license file /Users/zuck016/gurobi.lic
Set parameter TokenServer to value leghorn.emsl.pnl.gov
Initializing experiment deterministic_20210706.184055
/Users/zuck016/.pyenv/versions/anaconda3-2020.11/lib/python3.8/site-packages/cobra/util/solver.py:508: UserWarning: Solver status is 'infeasible'.
  warn(f"Solver status is '{status}'.", UserWarning)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-1-5f835dc253a1> in <module>
      8 #bsc_composite_det = BioscrapeCOBRAdeterministic(get_bioscrape_cobra_config())
      9 
---> 10 biocobra_data_det, comp0_det = simulate_bioscrape_cobra(
     11     initial_glucose=10, #mM
     12     initial_lactose=200, #mM

~/Projects/AgileBio/WholeCell/vivarium-notebooks/bioscrape_cobra/simulate.py in simulate_bioscrape_cobra(division, stochastic, spatial, initial_glucose, initial_lactose, initial_agent_states, bounds, n_bins, depth, diffusion_rate, jitter_force, divide_threshold, external_volume, n_agents, halt_threshold, total_time, sbml_file, emitter, output_type, parallel)
    599                 biocobra_experiment.update(sim_step)
    600     else:
--> 601         biocobra_experiment.update(total_time)
    602 
    603     # print runtime and finalize

~/.pyenv/versions/anaconda3-2020.11/lib/python3.8/site-packages/vivarium/core/engine.py in update(self, interval)
    644                     # TODO(chris): Is there any reason to generate a process's
    645                     #  schema dynamically like this?
--> 646                     update = self._process_update(
    647                         path, process, store, states, process_timestep)
    648 

~/.pyenv/versions/anaconda3-2020.11/lib/python3.8/site-packages/vivarium/core/engine.py in _process_update(self, path, process, store, states, interval)
    432         """
    433 
--> 434         update = self.invoke_process(
    435             process,
    436             path,

~/.pyenv/versions/anaconda3-2020.11/lib/python3.8/site-packages/vivarium/core/engine.py in invoke_process(self, process, path, interval, states)
    401             return self.parallel[path]
    402         # if not parallel, perform a normal invocation
--> 403         return self.invoke(process, interval, states)
    404 
    405     def _process_update(

~/.pyenv/versions/anaconda3-2020.11/lib/python3.8/site-packages/vivarium/core/engine.py in __init__(self, process, interval, states)
    188         self.interval = interval
    189         self.states = states
--> 190         self.update = invoke_process(
    191             self.process,
    192             self.interval,

~/.pyenv/versions/anaconda3-2020.11/lib/python3.8/site-packages/vivarium/core/engine.py in invoke_process(process, interval, states)
    124     """
    125     if process.update_condition(interval, states):
--> 126         return process.next_update(interval, states)
    127     return {}
    128 

~/.pyenv/versions/anaconda3-2020.11/lib/python3.8/site-packages/vivarium_cobra/processes/cobra_fba.py in next_update(self, timestep, states)
    270 
    271         # solve the fba problem
--> 272         objective_exchange = self.fba.optimize() * timestep  # mmol/L/s
    273         exchange_reactions = self.fba.read_exchange_reactions()
    274         exchange_fluxes = self.fba.read_exchange_fluxes()  # mmol/L/s

~/.pyenv/versions/anaconda3-2020.11/lib/python3.8/site-packages/vivarium_cobra/library/cobra_wrapper.py in optimize(self)
    339             self.solution = self.model.optimize()
    340 
--> 341         return self.objective_value()
    342 
    343     def external_reactions(self):

~/.pyenv/versions/anaconda3-2020.11/lib/python3.8/site-packages/vivarium_cobra/library/cobra_wrapper.py in objective_value(self)
    327     def objective_value(self):
    328         if self.solution:
--> 329             objective_value = self.solution.objective_value * self.flux_scaling
    330             return objective_value
    331         else:

TypeError: unsupported operand type(s) for *: 'NoneType' and 'float'

It seems that checking if self.solution should have caught the problem, but perhaps infeasible solutions are still solutions.

The simple fix would be to check if self.solution and self.solution.objective_value

vivarium-collective / vivarium-notebooks

Bioscrape and COBRA in a single deterministic cell results in infeasible solution if total_time > 360 #5