September 22, 2017, Friday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links


MBW:Oscillatory behavior in enzymatic control processes, rights and wrongs

From MathBio

Jump to: navigation, search

Summary

Starting with B.C. Goodwin's paper Oscillatory behavior in enzymatic control processes[1] from 1965 we investigate the two models presented in the paper using Simbiology in Matlab. A review of the paper is written by Scott Hoch. The general idea is that protein synthesis from mRNA is controlled via a negative feedback loop dependent on protein concentration. The paper presents two models simulating this, one for a single mRNA protein pair, and one where two pairs of mRNA and protein oscillators are coupled together. Using Simbiology we are able to reproduce the results from the paper, showing that the simulated negative feedback loop can create oscillations in both mRNA and protein concentrations. Also, when two oscillators are coupled together, we were able to produce beat frequencies in one of the oscillators, which could be regarded as a cells mechanism for long term concentration control, for example over a 24 hour run, where a certain protein can be down regulated during night and up regulated during day (or the other way around). We were also able to show that some of the results presented by Goodwin were incorrect, which also is reported in literature [3].

Introduction

Figure 1- mRNA creating a protein in the cell

Normal science is done by performing experiments (both physical and theoretical) within a paradigm, the set of assumptions that define a scientific discipline at any given time. These experiments are useful for discovering the details that are a consequence of the fundamental assumptions of the science. A more rare form of science is a paradigm shift occurs when a discovery forces the scientific community to change the framework from which they understand their field of science. In the 1960’s, many discoveries coupled with Brian Goodwin’s paper “Oscillatory behavior in enzymatic control processes” created a paradigm shift in the way that scientists understood the behavior of cells. Prior to the 1960’s, it was believed that molecules in the cell would move towards steady state concentrations. These concentrations determined by the environment inside and around the cell. Experimental data however was beginning to show evidence of biological “clocks” inside of living cells that would determine when processes would occur [2]. This evidence was an indicator that there must be oscillatory motion within a cellular organism and that scientists should change the way they thought about cells. One thing that was not understood and therefore hindering a paradigm shift was how a cell could keep track of many cyclic processes that have different cycle periods (sleep, photosynthesis, puberty, etc.) at the same time. In the 1965, Brian Goodwin’s paper mathematically described how two coupled biological oscillators could create a gamut of frequencies by varying the coupling constant. This provided a model for how two oscillators could control the timing for several biological processes all on different timescales. This discovery helped pave the way for a paradigm shift that matches our current understanding of both cellular organisms, and of coupled oscillators. This wiki will present a simplified model of a single biological oscillator and then couple two biological oscillators to show that different timescales can be created with such a system. The paper will be modeled after work done in Goodwin’s 1965 work.

Mathematical model

Single oscillator

For a single oscillator shown in Figure 2, mRNA produces proteins at a rate proportional to the number of mRNA in the cell. Therefore the only way to regulate the number of proteins in a cell, is to regulate the number of mRNA in a cell. For this reason the cell monitors the number of proteins it contains and then produces mRNA at a rate inversely proportional to the number of proteins in the cell. This is the negative feedback that will allow our system to oscillate. This behavior shows one mechanism by which a cell can keep track of time. This behavior is shown in the equations below where the variable X represents the concentration of mRNA, and Y represents the number of proteins in a cell. b and B are the death rates of mRNA and proteins.


{\begin{array}{c}(1)\ {\frac  {dX_{{1}}}{dt}}={\frac  {a_{{1}}}{A_{{1}}+k_{{1}}Y_{{1}}}}-b_{{1}}\\\\(2)\ {\frac  {dY_{{1}}}{dt}}=m_{{1}}X_{{1}}-B_{{1}}\end{array}}


with steady states given by

{\begin{array}{c}\\(3)\ X_{{1}}^{{*}}={\frac  {B_{{1}}}{a_{{1}}}}\\\\(4)\ Y_{{1}}^{{*}}={\frac  {(\left(a_{{1}}/b_{{1}}\right)-A_{{1}})}{k_{{1}}}}.\end{array}}

Coupled oscillators

There are many types of mRNAs all making different types of proteins simultaneously in a living cell. With so many proteins being created by and inhibiting the production of mRNA, it is to be expected that these oscillatory cycles should interact with each other. For example a high population of a protein can inhibit the production of an mRNA that is designed to make a different protein. This behavior can be modeled as two oscillators each creating negative feedback for itself and for the other oscillator. Such a system is known as a coupled oscillator. Though there are many proteins mRNA pairs in any given cell, we begin by examining the simplest possible system of two coupled oscillators. The behavior of these two oscillators is given by


{\begin{array}{c}(5)\ {\frac  {dX_{{1}}}{dt}}={\frac  {a_{{1}}}{A_{{1}}+k_{{11}}Y_{{1}}+k_{{12}}Y_{{2}}}}-b_{{1}},{\frac  {dY_{{1}}}{dt}}=m_{{1}}X_{{1}}-B_{{1}}\\\\(6)\ {\frac  {dX_{{2}}}{dt}}={\frac  {a_{{2}}}{A_{{2}}+k_{{22}}Y_{{2}}+k_{{21}}Y_{{1}}}}-b_{{2}},{\frac  {dY_{{2}}}{dt}}=m_{{2}}X_{{2}}-B_{{2}}\\\end{array}}

where k12, and k21 describe the coupling of the two oscillators. If these two variables are zero, then the oscillators are no longer coupled.

Results and analysis

Single oscillator

The model

We reproduced the model for the simple oscillator of mRNA and protein presented in the paper using Simbiology in Matlab. The model can be seen in figure XX. As we can see, the mRNA compartment and the protein compartment are not directly coupled together. This is because the reaction set is not of a mass action type where two reactants come together to form a product. Rather, the synthesis of mRNA is inversely proportional to the concentration of protein, and for the purpose of this model, the protein is solely produced by the mRNA independent of protein concentration. This is what the model simulates.

Figure 2- Simbiology model of a single negative feed back oscillator.

Initially we set all parameters as presented in the paper. These were


a=72,A=36,k_{{1}}=1,b_{{1}}=2,m_{{1}}=1,B_{{1}}=0,X_{{1}}(t=0)=7,Y_{{1}}(t=0)=10.

The only change from Goodwin’s parameters and initial conditions was to the initial amount of protein, which in the paper was set to -10 for scaling purposes. This does clearly not make sense in biology, so we chose to set the initial protein concentration to 10 instead. The results from running the model is shown in Figure 3. As expected we see that both species

Figure 3- Nonlinear oscillations of a single biological oscillator.
Figure 4- Oscillations are created from initial concentrations near the steady state (0,0).

oscillate with the protein amount lagging slightly behind the mRNA amount. This makes sense because in order to increase the number of protein in the cell, the mRNA concentration must rise first. Also, as reported in the paper, the oscillations of the protein appear sinusoidal, whereas the oscillations of the mRNA are nonlinear. The steady state of this model with parameters as defined above has solutions of both mRNA and protein concentrations equal to zero. We tried to investigate this steady state by running a simulation with initial mRNA amount at 0.1 and initial protein amount at 0. As seen in Figure 4, the concentrations moved away from the steady state, rather than towards it as one might expect. These oscillations occurred with a smaller amplitude than shown in Figure 3. Also, the mRNA in this simulation came out as nearly sinusoidal (linear) oscillations meaning the system is not being over driven. The fact that oscillations are large even though the initial conditions are very close to the steady state variables, indicates that the steady state is only weakly stable.

Phase plot

We used pplane (program which allows you to write in two ODEs and plots out the phase plane. It is free and is easiest found through a simple Google search) to produce phase plots of protein versus mRNA of the first simulation described above. This gave orbits around (0,0) as expected. Again, phase plots very near to the steady state solution gave orbit solutions and not spirals towards (0,0), which also indicates that to reach steady state, neither mRNA nor protein can be present at all, and that it is a weak stable steady. The phase plot is shown in figure 5.

Figure 5- Phase plot of mRNA versus protein for a single oscillator.

Coupled oscillators

Figure 6- Schematic diagram and simbiology model of two coupled oscillators

Now that we are convinced we can model a single oscillator successfully, we present the coupled oscillator model using Simbiology. The model can be seen in figure 6. Again, none of the species in the model (mRNA1, mRNA2, protein1 and protein2) are coupled together in the diagram model, since these coupled oscillators do not follow simple mass action dynamics either, as discussed for the single oscillator. To run the simulations we used the same parameter values as presented in the paper. These were


{\begin{array}{l}\\a_{{1}}=360;A_{{1}}=36;b_{{1}}=10;k_{{11}}=1;k_{{12}}=0;m_{{1}}=0.5;B_{{1}}=0;\\\\a_{{2}}=360;A_{{2}}=43;b_{{2}}=10;k_{{22}}=1;k_{{21}}=0;m_{{2}}=0.6;B_{{2}}=0\\\\\end{array}}

Figure 7- Two uncoupled oscillators with different resonant frequencies

With these parameters the steady states of both oscillators are (0,0), as in the single oscillator case. Since both coupling constants, k12 and k21, are set to zero, the system should also now act as two uncoupled oscillators. Also note that the only difference between the two oscillator’s parameters is their values of m. The first system has m1 = 0.5, and the second system has m2 = 0.6. From equations 5 and 6 we see that a higher value of m would give a higher rate of change in protein, and therefore we can expect oscillations with higher frequency in oscillator pair 2 than oscillator pair 1. Running a simulation on the model confirms this, as can be seen in Figure 7.

Figure 8- The faster oscillator driving the slower oscillator creates a beat frequency in the slow oscillator.

We now couple the oscillators together by giving k12 and k21 values different than zero. All other parameter values are unchanged. We first set k12 to 0.3 and let k21 stay at zero. This means that oscillator two impacts oscillator one via negative feedback, while oscillator one has no influence on oscillator two (the faster oscillator can impact the slower). The result is shown in figure 8. As expected, no change is observed in oscillator two, it has the same amplitude an frequency as before. Beats are on the other hand observed in oscillator one. The frequency within the beats is the same as for the uncoupled oscillators


Figure 9- The slower oscillator driving the faster oscillator creates a beat frequency in fast oscillator.
Figure 10- Both oscillators exert negative feedback on each other.
Figure 11- This coupled oscillators do not have a stable orbit about the point (0,0).

Reversing the coupling such that oscillator one can impact oscillator two but not the other way around (with k12=0 and k21=0.3) is shown in Figure 9. Now no change is seen in oscillator one, but beats are observed in oscillator two. The beat frequency is much lower in this case than when the coupling was the other way around. In this last simulation, the slower oscillator is allowed to impact the faster one. It therefore makes sense that the beat frequency, which is a result of oscillator one impacting oscillator two, is lower in this case than when the faster oscillator drives the slower. It is interesting to note here that the beat frequency in Figure 9 is ~22 times faster than the oscillation frequency of the oscillator itself. This leads to a mechanism for keeping track of both hourly and daily cycles (24 hours per day) The next problem we address is what happens if both oscillators are coupled such that they both exert negative feedback on themselves, and on the other oscillator. We begin by setting k12 = 1.2 and k21 = 0.77. This means that oscillator one (the slower one) impacts oscillator two stronger than the other way around. Figure 10 shows the result. We see that protein levels of both pairs oscillate with the same frequency and approximately the same amplitude. Also, we see that neither of the plots are sinusoidal, but are highly impacted by the original frequencies of the uncoupled oscillators. As we did for the single oscillator, we now study the stability of the steady state of the coupled oscillators. We do this by setting the initial value of protein 1 to 0.1, and the initial values of the three other species to zero. All other parameter values were set as in the previous simulation, including values for k12 and k21. If the steady state is highly stable, we would expect the solution to drop to zero and stay there. As shown in Figure 11, this is not what happens. Instead of dropping to zero, all four species immediately go into stable, non-zero oscillations. As before, this show us that the steady state at (0,0) for both oscillators are weakly stable.

Discussion

Negative concentrations of species

The first issue which must be addressed is the scaling of the concentrations in the figures shown above. A lot of the concentrations shown have negative values while that is not physically possible. This is a result of the parameter values chosen by Goodwin to allow the concentrations to stay on the screen scale of the computer he was using for numerical simulations. This is not a concern to us because the purpose of this paper is to show that oscillatory behavior in cells can exibit multiple frequencies when there are more than one oscillators coupled to each other. The figures shown above succeed in doing this.

Agreement with Goodwin

The models which were created in Simbiology agree with the findings of Goodwin's 1965 paper in all cases except for the case of the faster oscillator exerting negative feedback on the slower oscillator. Our study found that this coupling scheme created beat frequencies in the slower oscillator and not in the faster oscillator. This disagreed with Goodwin's simulation that showed the coupling caused the two oscillators to sync to the same frequency. Upon farther investigation, we found that the findings of Goodwin were incorrect due to errors in programming by the people at MIT who ran the simulations for him [3]. To see Goodwin's results, please see the wiki page on his paper[1].

Questions for biologists

This paper presented several models that could describe the way that mRNA protein oscillators could be coupled together, but gives no information on which coupling mechanism is more probable. This question could be answered by an experiment in the following way. It is seen that when there is mutual coupling between the oscillators that the total amount of protein in the cell is conserved, meaning that the oscillations of each protein are out of phase. However, when only one protein is allowed to exert negative feedback on the other, the driving oscillator has a constant amplitude while the driven oscillator has a beating amplitude. A measurement of the total concentration of proteins in a cell in time could definitively determine which of these models describes nature in a more accurate way.

Conclusion

In conclusion, we have shown that for the most part, the simulations preformed by Goodwin in the 1960's were preformed correctly and showed how coupled oscillators can create a gamut of frequencies that can be used for biological clocks. Furthermore, we were able to take the models proposed by Goodwin and came up with an experiment that would be able to show which model exhibits behavior that most closely replicates nature. We also agreed on the fact that having zero concentrations of proteins and mRNA is an unstable equilibrium. This is a good thing because if the cell pushed concentrations of proteins to zero, the cell would no longer be able to exist.

References

1. B. C. Goodwin. Oscillatory behavior in enzymatic control processes. Adv. Enzyme Regula tion, (1965) 3:425�-4384 438.


2. J. L. CLOUDSLEY, ON, Rhythmic Activity in Animal Physiology and Behaviour, Academic Press, New York and London (1961).


3. J. S. GRIFFITH. Mathematics of Cellular Control Processes. I. Negative Feedback to One Gene. J. Theoret. Biol. (1968) 20: 202-208.

Matlab code

This Matlab code were used in the simulations. It was created using Matlab's Simbiology.

Single oscillator

function data = runsimulation(modelobj)

% Get the active configuration set.
cs = getconfigset(modelobj, 'default');

% Run simulation.
data = sbiosimulate(modelobj, cs, []);

% Plot the results.
plottype_Time(data, '<all>', 'one axes');

% ----------------------------------------------------------
function plottype_Time(tobj, y, plotStyle)
%TIME Plots states versus time.
%
%    TIME(TOBJ, Y, PLOTSTYLE) plots the results of the simulation for the
%    species with the specified Y versus time. If PLOTSTYLE is 'one axes'
%    then data from each run is plotted into one axes. If PLOTSTYLE is
%    'subplot' or 'trellis' then data from each run is plotted into its own
%    subplot. If Y is '<all>' then all data will be plotted.
%
%    See also GETDATA, SELECTBYNAME.

if (length(tobj) > 1)
    switch (plotStyle)
    case 'one axes'
        sbioplot(tobj, @timeplotdata, [], y);
    case 'subplot'
        % Include a legend with the subplot.
        sbiosubplot(tobj, @timesubplotdata, [], y, true);  
    case 'trellis'
        sbiotrellis(tobj,@timesubplotdata, [],y);   
    end
else
    % Plot Data.
    timesubplotdata(tobj, [], y);
    
    % Label the plot. 
    xlabel('Time');
    ylabel('States');
    title('States versus Time');
    
    % If Y is '<all>' get all the data names.
    if strcmpi(y, '<all>')
        [~, ~, names] = getdata(tobj);
    else
        names = y;
    end
    
    % Create legend.
    leg = legend(names, 'Location', 'NorthEastOutside');
    set(leg, 'Interpreter', 'none');
end

%--------------------------------------------------------
function [handles, names] = timeplotdata(tobj, ~, y)

colors    = get(gca, 'ColorOrder');
numColors = length(colors);

% Preallocate handles
if strcmpi(y, '<all>')
    [~, ~, names] = getdata(tobj(1));
else
    [~, ~, names] = selectbyname(tobj(1), y);
end
handles = zeros(length(names), length(tobj));

for i=1:length(tobj)
    % Get the data from the next run.
    nexttobj = tobj(i);

    % Get the data associated with y.
    if strcmpi(y, '<all>')
        [time, data, names] = getdata(nexttobj);
    else
        [time, data, names] = selectbyname(nexttobj, y);
    end

    % Error checking.
    if size(data,2) == 0
        error('Data specified do not exist.');
    end
    
    % Plot data. If there is only one state use different colors for runs.
    if(size(data,2) ==1)
        hLine = plot(time, data, 'color',colors(mod(i-1,numColors)+1,:));
    else
        hLine = plot(time, data);
    end
    handles(:,i) = hLine;
end

% ---------------------------------------------------------
function timesubplotdata(tobj, ~, y)

% Get Data to be plotted.
if strcmpi(y, '<all>')                
    [time, data] = getdata(tobj);
else
    [time, data] = selectbyname(tobj, y);
end

% Error checking.
if size(data,2) == 0
    error('Species specified do not exist.');
end

% Plot Data.
plot(time, data);

Coupled oscillators

function data = runsimulation(modelobj)

% Get the active configuration set.
cs = getconfigset(modelobj, 'default');

% Run simulation.
data = sbiosimulate(modelobj, cs, []);

% Plot the results.
plottype_Time(data, {'unnamed.X1','unnamed.X2'}, 'one axes');

% ----------------------------------------------------------
function plottype_Time(tobj, y, plotStyle)
%TIME Plots states versus time.
%
%    TIME(TOBJ, Y, PLOTSTYLE) plots the results of the simulation for the
%    species with the specified Y versus time. If PLOTSTYLE is 'one axes'
%    then data from each run is plotted into one axes. If PLOTSTYLE is
%    'subplot' or 'trellis' then data from each run is plotted into its own
%    subplot. If Y is '<all>' then all data will be plotted.
%
%    See also GETDATA, SELECTBYNAME.

if (length(tobj) > 1)
    switch (plotStyle)
    case 'one axes'
        sbioplot(tobj, @timeplotdata, [], y);
    case 'subplot'
        % Include a legend with the subplot.
        sbiosubplot(tobj, @timesubplotdata, [], y, true);  
    case 'trellis'
        sbiotrellis(tobj,@timesubplotdata, [],y);   
    end
else
    % Plot Data.
    timesubplotdata(tobj, [], y);
    
    % Label the plot. 
    xlabel('Time');
    ylabel('States');
    title('States versus Time');
    
    % If Y is '<all>' get all the data names.
    if strcmpi(y, '<all>')
        [~, ~, names] = getdata(tobj);
    else
        names = y;
    end
    
    % Create legend.
    leg = legend(names, 'Location', 'NorthEastOutside');
    set(leg, 'Interpreter', 'none');
end

%--------------------------------------------------------
function [handles, names] = timeplotdata(tobj, ~, y)

colors    = get(gca, 'ColorOrder');
numColors = length(colors);

% Preallocate handles
if strcmpi(y, '<all>')
    [~, ~, names] = getdata(tobj(1));
else
    [~, ~, names] = selectbyname(tobj(1), y);
end
handles = zeros(length(names), length(tobj));

for i=1:length(tobj)
    % Get the data from the next run.
    nexttobj = tobj(i);

    % Get the data associated with y.
    if strcmpi(y, '<all>')
        [time, data, names] = getdata(nexttobj);
    else
        [time, data, names] = selectbyname(nexttobj, y);
    end

    % Error checking.
    if size(data,2) == 0
        error('Data specified do not exist.');
    end
    
    % Plot data. If there is only one state use different colors for runs.
    if(size(data,2) ==1)
        hLine = plot(time, data, 'color',colors(mod(i-1,numColors)+1,:));
    else
        hLine = plot(time, data);
    end
    handles(:,i) = hLine;
end

% ---------------------------------------------------------
function timesubplotdata(tobj, ~, y)

% Get Data to be plotted.
if strcmpi(y, '<all>')                
    [time, data] = getdata(tobj);
else
    [time, data] = selectbyname(tobj, y);
end

% Error checking.
if size(data,2) == 0
    error('Species specified do not exist.');
end

% Plot Data.
plot(time, data);