Guided example: lac operon as a difference equation model

Let’s revisit our example of the lac operon. For didactic reasons, we will keep the model components the same as in the Boolean model description of the system, i.e. all model components represent the activity of the biological entity modelled. However, we now describe their behaviour using quantitative variables that evolve over discrete time steps based on kinetic assumptions. We will model the activating and inhibiting interactions between the components with abstracted rate laws.

But remember, this is by far not the only option to model this system. We could also formulate the system using explicit biochemical reactions and derive equations via mass action kinetics (e.g. modelling the binding of cAMP to CAP (Catabolite Activator Protein) as: CAP_inactive + cAMP ⇌ CAP_active), which would require additional components.

To keep the model both simple and illustrative, we make the following modelling assumptions:

• All model components have a constant production and linear degradation rate, unless production or degradation is regulated by another component.
  
• Glucose and allolactose are treated as external inputs with constant levels.
  

Fixed model inputs

We assume that external glucose and internal allolactose concentrations are fixed in our scenario. They can be considered constants, though formally we could write the corresponding difference equations:

glucose (t + Δt) = glucose (t)

allolactose (t + Δt) = allolactose (t)

Dynamic model components

For the other model components we need to determine the equations f(x) that describe by how much each model component changes per unit of time.

cAMP: Glucose inhibits cAMP production and therefore, cAMP levels should increase when glucose level is low. If we don’t know any details about the mechanism of such a regulation, a common assumption is that the unknown process has a maximum rate (cAMP/unit time) and this rate is modulated by the glucose level in an inversely proportional way. 

For the modelling consumption to be linear and proportional to the current cAMP level we get:

CAP: The activity of CAP increases as cAMP levels rise because cAMP binding activates CAP. We assume there is a maximum rate at which active CAP can be produced, denoted by k_CAP,max. The actual production rate scales between zero and this maximum, depending on the concentration of cAMP. Specifically, when cAMP is absent, the production rate is zero, and when cAMP levels are very high, the rate approaches its maximum value. This maximum reflects factors such as the total amount of inactive CAP available (which is not explicitly modelled here) and the binding free energy of cAMP to CAP. A typical function capturing this behaviour is :

This function equals 0 when cAMP is 0, approaches 1 as cAMP becomes large, and reaches half its maximum value when cAMP equals K_CAP.

Again assuming the consumption term to be linear we get:

lac repressor: The lac repressor activity is inhibited by allolactose. We assume the active repressor protein is produced with a constant rate. Its consumption rate however is modulated by the allolactose level to increase with increasing levels of allolactose. The consumption rate should also be proportional to the active repressor level. A function that fulfills these requirements is:

This leads to: 

Polycistronic lac mRNA: Transcription of the lac operon is only possible when the repressor is not bound. Transcription is enhanced by binding of active CAP. Therefore, the production rate of lac mRNA decreases with increasing repressor levels and increases with increasing CAP levels. We can model the production rate using the following function:

Whereas the consumption rate is modelled proportional to the level of mRNA.

You might wonder: how can the production rate be represented as a continuous value when, in reality, the repressor is either bound or unbound, and the same goes for CAP? This would suggest that transcription is simply on or off.

The explanation is that although individual binding events are indeed binary, when we measure outputs like mRNA levels or fluorescence, we are observing averages across large cell populations or over time. These averages smooth out the on – off nature of single events, resulting in effectively continuous production rates.

Lac enzymes: We assume lac enzymes are translated proportionally to mRNA levels and degraded proportionally to protein concentration:

Lactose metabolism: The model component lactose metabolism is a proxy for a downstream metabolic activity. Like the readout of a spectrophotometric enzyme assay that measures fluorescence absorbance of a metabolic reaction product but not the enzyme levels directly. A simple way to represent this is:

lacMetabolism = k_conv​⋅lacEnzymes

The conversion factor k_conv could be determined experimentally (e.g. via a calibration curve in a fluorescence assay).

Using the established difference equations, the model updates each component at every time step by first calculating the change based on the current values, then adding this change to the current values to obtain the new state. These updated values are then used to determine the changes in the following time step.

Exploring model behaviour

To explore the model dynamics we simulate the model by computing step by step updates of all model components according to their difference equations. Figure 14 shows the model dynamics for a scenario in which glucose is absent and allolactose is present over a duration of 10 time steps.