PID Control Demonstration

Those of you who are interested in feedback control may have come across the infamous PID controller. To many engineers, PID control is the answer to all their problems. What is it for? PID controllers are a sophisticated way to stabilize a system that is under continual threat of disturbance. The simplest example is the cruise control in a car. Most of us probably set the cruise control and forget about it. What’s happening behind the scenes is however much more interesting. As we go up a hill we notice the engine revving up to maintain the set cruise speed and when we go down a hill the engine slows down. The simplest feedback mechanism is where the controller compares the speed of the car with the set cruise speed and adjusts the engine accordingly. The secret to a good controller such as this is what it does with the error it finds between the actual speed of the car and the desired speed. The simplest solution is to adjust the speed of the engine directly in proportion to the difference, this is called proportional control and is most familiar to engineers. The problem with proportional control is twofold:

1. Too much proportional control and the response overshoots, perhaps so much so that it goes into uncontrolled oscillations, not a good thing in a car.

2. You usually have to set the cruise control slightly higher than the speed you actually want because a proportional controller finds it difficult to actually match the set speed unless the proportional control is strong enough (See problem 1!)

Both these problems can be solved by combining proportional control (The P in PID) with two others, integral (The I in PID) and derivative control (The D in PID). Derivative control is perhaps the easiest to understand. What derivative control does is measure how fast the error between desired and actual engine speed changes. If the difference is changing too fast then we want to slow down the response to avoid overshooting. Derivative control will therefore help get rid of the instability (Too much derivative control can however also result in instability).

The second problem is that the speed of the car doesn’t actually settle to the desired speed. This is fixed by using integral control. This is where the error is continually added up to form a net error and it is the net error than is then used to control the speed of the engine. So long as there is even a little bit of error (ie Desired speed 50 mph, actual speed 49 mph), the net error will increase and drive the car speed to the desired speed until the actual error is zero. Integral is particularly useful if there is a continual disturbance at work, for example, a headwind.

To demonstrate these capabilities, I have written a small PID control demonstration tool (for windows only I am afraid) that simulates each of the three elements in a PID controller. The screenshot below shows the application. A user can control the strength of the different controlling elements or remove them all together. To test the response a user can apply a step change to the desired car speed and see how the system responds. The simulation is continuous so it is easy to see the effects.

Download the PID Controller Demo here

Partly supported by the NSF

Front-Loading of Flux Control

Originally Posted on December 20, 2019 by hsauro

There is an interesting property that straight chain pathways have which I call front-loading. The phenomenon has been known for some time and I describe it in detail in my recent textbook on Metabolic Control Analysis. Let’s say we have a straight chain of reactions, something like:

  

There are no long-range feedback loops (other than short-range product inhibition) and it is assumed that for a given reaction, increases in substrate cause the reaction rate to increase and increases in the product causes the reaction to decrease. Each reaction is reversible using a simple reversible mass-action rate law: . I assume that  and  are fixed species. Given suitable rate constants and equilibrium connstants for each reaction, the system will admit a steady state. A question to ask is what is the distribution of flux control across the chain? There are two ways to answer this, we can use theory, this gives us what the distribution is and insight into why or we can do simulation but this will just tell us what the distribution is likely to be. For now let’s do a simulation.

The code is shown below. The code defines a four-step straight-chain, runs 10,000 simulations, randomizing the rate constants for each simulation. Each time it does a simulation we compute the flux control at each step, store this information and then at the end, plot histograms of the distribution of the flux control coefficients. Scroll down to see the plots. The average values for the flux control coefficients is shown in the histogram:

The distribution of flux control is shown below. The vertical axis is the frequency and the x-axis the value of the flux control coefficient from zero to 1.0. The red curve corresponds to the flux control for the first step, on average this has the highest control. The yellow distribution corresponds to the last step of the pathway, you can see it is least likely to have any significant flux control. The conclusion is that given a straight chain with a random set of parameters, on average the first step will have the highest control, progressively decreasing as we work our way down the pathway. Why is this? The explanation is a thermodynamic one which in turn boils down to have easy it is for perturbation to travel up and down the chain. So long as the thermodynamic gradient is from left to right, it is easier for perturbations to propagate downstream compared to upstream. Since flux control is really a measure of how much a perturbation has on the steady-state flux, the easier a perturbation can travel the higher the flux control.

This is not to say that it isn’t possible to get excellent flux control in downstream steps, it’s just that fewer parameter combinations will achieve that.

# Monte Carlo simulation of a straight chain pathway
# Samples parameter values while keeping Keq constant
# Plots the distribution of control coefficients
# Not the most elegant code but I wrote it quickly
 
import tellurium as te
import roadrunner
import random
import pylab as pl
 
r = te.loada("""
   J1: $Xo -> S1; k1*(Xo - S1/Keq1);
       S1 -> S2;  k2*(S1 - S2/Keq2);
       S2 -> S3;  k3*(S2 - S3/Keq3);
       S3 -> $X1; k4*(S3 - X1/Keq4);
     
     k1 = 0.1; k2 = 0.1;
     k3 = 0.1; k4 = 0.1;
     
     Keq1 = 4;
     Keq2 = 3;
     Keq3 = 2;
     Keq4 = 1;
     Xo = 5;
     X1 = 0.1;
""")
 
m = r.simulate (0, 10, 100);
 
aC1 = 0; aC2 = 0; aC3 = 0; aC4 = 0;
aC1a = []; aC2a = []; aC3a = []; aC4a = [];
n = 1000
a = 0
upperLimitK = 50
for i in range (0,n):
    r.setValue ('k1', random.uniform(0, upperLimitK))
    r.setValue ('k2', random.uniform(0, upperLimitK))
    r.setValue ('k3', random.uniform(0, upperLimitK))
    r.setValue ('k4', random.uniform(0, upperLimitK))
    
    try:
      r.simulate()    
      r.steadyState()
      C1 = r.getCC ('J1', 'k1')
      C2 = r.getCC ('J1', 'k2')
      C3 = r.getCC ('J1', 'k3')
      C4 = r.getCC ('J1', 'k4')
      aC1 = aC1 + C1
      aC2 = aC2 + C2
      aC3 = aC3 + C3
      aC4 = aC4 + C4
      aC1a.append (C1)
      aC2a.append (C2)
      aC3a.append (C3)
      aC4a.append (C4)
    except:
        a = a + 1
 
print (aC1/n, aC2/n, aC3/n, aC4/n) 
 
bins = 100
 
pl.hist(aC1a, bins=bins, histtype='stepfilled', density=True, color='r', alpha=0.5, label='C1')
pl.hist(aC2a, bins=bins, histtype='stepfilled', density=True, color='b', alpha=0.5, label='C2')
pl.hist(aC3a, bins=bins, histtype='stepfilled', density=True, color='g', alpha=0.5, label='C3')
pl.hist(aC4a, bins=bins, histtype='stepfilled', density=True, color='y', alpha=0.5, label='C4')
print ("fails = ", a)

Stream Plots Using Tellurium

Originally posted on December 18, 2019 by hsauro

Every wanted to do a phase plot for a toggle switch or any two-dimensional model? Here is an example of using a stream plot combined with SBML models using Tellurium based on a simple toggle switch.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46

import numpy as np import matplotlib.pyplot as plt import tellurium as te   Y, X = np.mgrid[0:8:200j, 0:8:200j]   r = te.loada(''' J1: -> x; k1/(1+y^n1) - k2*x; J2: -> y; k3/(1+x^n2) - k4*y;   x = 4; y = -4; k1 = 12; k3 = 12; k2 = 2; k4 = 2 n1 = 4; n2 = 4 ''')   m = r.simulate (0, 8, 100)   U = np.empty_like(X) V = np.empty_like(X)   for (i,j), value in np.ndenumerate(X):     r.x = X[i][j]; r.y = Y[i][j]     U[i][j] = r.J1; V[i][j] = r.J2   plt.figure(figsize=(10,7)) plt.xlabel('x', fontsize='16') plt.ylabel('y', fontsize='16') plt.ylim((-0.1,8)) plt.xlim((-0.1,8)) plt.streamplot(X, Y, U, V, density=[2, 2])   r.steadyStateSolver.allow_presimulation = False   r.x = 5.999; r.y = 0.00462 r.steadyState() plt.plot(r.x, r.y, 'ro', markersize=12)   r.x = 0.0046; r.y = 5.999 r.steadyState() plt.plot(r.x, r.y, 'bo', markersize=12)   r.x = 1.359; r.y = 1.359 r.steadyState() plt.plot(r.x, r.y, 'go', markersize=12)   plt.show()

The colored circles mark the three steady states. The green circle in the middle is on a saddle-node which is unstable (unless you’re exactly on the saddle ridge) and the two on either side are stable nodes.