`ThinLayerPysces`: proteins not updated during simulation

[jmrohwer]

I have been testing ThinLayerPysces with a complex dataset comprising 15 measurements and 2 cascade reactions. The parameters I have obtained with PyEnzyme differ from those I obtained when coding the fit by hand using lmfit. It turns out that the catalyst (protein) concentrations differ between measurements. These are specified in the rate laws of the two reactions (as modifiers). But it appears that the protein values are not updated for the various measurements, see here:

PyEnzyme/pyenzyme/thinlayers/psyces.py

Lines 530 to 551 in 8a29bb7

	def to_pysces_model(self, model: pysces.model):
	"""
	Apply initial conditions to a PySCeS model.
	This method sets the simulation time and initial species concentrations
	in the PySCeS model, ensuring zero values are replaced with small positive values.
	Args:
	model (pysces.model): The PySCeS model to update.
	Returns:
	pysces.model: The updated model with initial conditions set.
	"""
	model = dill.loads(dill.dumps(model))
	model.sim_time = np.array(self.time)
	model.__dict__.update(
	{
	f"{species_id}_init": value if value != 0 else 1.0e-9
	for species_id, value in self.species.items()
	}
	)
	return model

Only the initial values of the species in the measurements are updated. I guess the values of the constant species (if there are any), as well as of the proteins, should be updated as well as these may differ between measurements. In the PySCeS PSC files these enter as "fixed species".

[jmrohwer]

@JR-1991 any update on this? We discussed this briefly at the EnzymeML workshop last week. I need the functionality to work on the dataset for Torsten's NMR paper. EnzymeML document used is attached for reference.

[jmrohwer]

pgm-eno.json

[JR-1991]

@jmrohwer, I’ve created a PR that addresses the issue. Here’s a screenshot of the initial conditions and measurements for the PGM-ENO dataset. The protein initial is now included in the ICs:

[jmrohwer]

@JR-1991 unfortunately the code in #90 still does not give the correct answer. The issue is with the code quoted above. Specifically:

model.__dict__.update( 
             { 
                 f"{species_id}_init": value if value != 0 else 1.0e-9 
                 for species_id, value in self.species.items() 
             } 

The problem is that PySCeS model species only have an _init attribute if the species is part of the stoichiometry (and thus changes during the simulation), i.e. it is a variable in the simulation. Modifiers (e.g. enzyme concentrations, protein concentrations), if they are not part of the stoichiometry, are treated as "parameters" from a modelling point of view. However, they might change from measurement to measurement.

The consequence is that the above code will write model.prot_init but that does absolutely nothing in PySCeS as model.prot is not a variable species. So the model object gets a new attribute (model.prot_init) that does nothing and just sits there. In this case model.prot should rather have been updated itself.

The way to handle this is:

• filter InitMap.species as before to ensure that the species remaining occur in a rate equation somewhere
• for each of these species:
  • if it is a variable species (in the PySCeS sense of the word), then update model.<species_id>_init as above. This can easily be checked because all variable model species are in a list model.species which is an attribute of the PySCeS model.
  • if the species is not part of the variables, then update model.<species_id> directly. This would be the case for modifiers, enzyme concentrations, protein concentrations etc.

I hope the explanation is clear...

[jmrohwer]

In addition, I would change

f"{species_id}_init": value if value != 0 else 1.0e-9 

to

f"{species_id}_init": value 

We should not be second-guessing the intentions of the modeller. If the modeller wants to replace zero with a small value to avoid numerical issues, it would be upon them to do this themselves.

For example, consider a model with concentrations of nM (realistic for hormones, for example). However, the model itself could have been specified in M. Then a value of 1e-9 is a perfectly acceptable value for an initial concentration.

[JR-1991]

@jmrohwer thanks for clarification. Just to make sure I understand this correctly, in order to have this on a per-measurement basis, the protein concentration should be passed without the _init affix? If so, the fix would look something like the following?

for species, value in self.species.items():
      if hasattr(model, f"{species}_init"):
          setattr(model, f"{species}_init", value)
      elif hasattr(model, species):
          setattr(model, species, value)
      else:
          raise ValueError(f"Species {species} not found in model")

I have removed the 1e-9 assumption here and would also commit it. Makes total sense to not assume things.

[JR-1991]

Using this fix resulted in the following fit. Is this expected?

[jmrohwer]

@JR-1991 This indeed looks correct but I cannot reproduce it on my side and the changes don't seem to be present in the code. Did you push the changes to #90 ? If not, please do.