Hello again.

We are now at the third part

of our lecture on Introduction to Dynamical Systems.

In the last lecture in part two of the Dynamical Systems Lecture Series, we

discussed how you would numerically solve systems

of ODE's specifically we discussed using Euler's method.

Now in this lecture we're going to talk about not just

how you solve them numerically but how you analyze these ODE

systems to get, to understand general properties of their system, and

we're going to focus on this concept of stability of ODE systems.

So first we'll discuss what we mean in general when we talk about stability.

Then I'm going to show a couple of one-dimensional examples

so you can understand how, how stability works in that context.

And then we're going to introduce

phase-plane techniques for analyzing two-dimensional systems.

Although most of what we explained about

these phase-plane techniques are going to be

covered in more depth in the, in the fourth part of this series of lectures.

And then I'm going to briefly introduce an example that we're going to deal

with in in part four which is

a mathematical model of glycolytic oscillations in yeast.

The focus of this lecture as mentioned is

going to be an analyzing stability of OED systems.

So to introduce what we mean when we talk about stability, let's go

back once again to this generic three

component repressive network that we've seen before.

Where you have three protein species, a, b, and c,

that could either be dephosphorylated or they could be phosphorylated.

And a, b, and c regulate either the

dephosphorylation or the phosphorylation reactions of the other species.

And as we have discussed, this scheme implies a set of differential

equations, which we've gone through and looked at as follows over here.

What we haven't discussed as much are the model parameters in the system.

And when we look at this we can see several constants in here.

That are that are model parameters that can be varied.

And that can govern the behavior of the system.

One is we have Kcat and km's for the phosphorylation reactions.

And those are the ones that are over here for, labeled as kinases.

And Kcats are the little k's in the numerator, and then the

km's are the big k's for

the phosphorylation reactions in the denominator.

Similarly we have Kcats and Kms for the dephosphorylation reactions, those

are the analogous parameters over here on these left, these leftmost terms.

And then the other things we have, the other parameters we have our total amounts

of A B and C that is A total B total and C total here.

And as we have discussed, these

equations are solved using standard numerical techniques.

And I'll just give a couple examples of what you

can see when you vary the parameters in this case.

There are two two general categories of solutions to this system.

With one set of parameters, we can see that A

increases to a steady level, B goes to a different steady

level that's much lower than the steady level of A,

and then C goes to a level that's very, very low.

But if you continue to run this for longer and longer periods

of time, it would continue to stay at this very steady level.

Once A gets to this level here, right around

0.9, it stays there, and it will stay there forever.

But with a different parameter set, we see this sort of behavior.

We can see the sustained oscillations able to go

up and then go down and similarly B and

C will also increase and decrease as a function

of time and these oscillations will go, go on forever.

So what we conclude from this in general is

that the parameter values can greatly influence the system behavior.

We have qualitatively different behavior over here on the left.

Where A, B, and C go to steady levels.

Compared to on the right, where A, B, and C continue to oscillate.

How do we understand these different behaviors?

Well, that's where the tools of dynamical systems come in.

And that's how we can, we can understand and

categorize the different types of behaviors that we observe.

[SOUND]

To illustrate how the tools of dynamical systems can be

used to analyze stability, we will consider a one-dimensional example.

This is a standard mathematical model of an

isolated cardiac myocyte, where different ionic currents are responsible

for changing the, the cells transmembrane potential, the

details of these ionic currents are not important here.

What is important here is that the

differential equation describing voltage across the cell membrane.

DV, dt is a negative of the ionic current.

That means the sum of, of the current through all of

these different channels, and pumps and transporters, divided by the capacitance.

And we want to do a, a relatively

simple experiment with this isolated cardiac myocyte model.

We begin with the myocyte at rest, which is approximately minus 85 millivolts.

Then we instantaneously change the voltage.

When we instantaneously change the voltage we

calculate the instantaneous ionic current, I ion.

And if we, once we calculate the instantaneous ionic current, then

we know what the instantaneous change is in the in the derivative.

So, we can fit, we can vary voltage to

whatever we want it to be instantaneously from rest.

And then we want to calculate what's the resulting

derivative going to be, what's the dV/dt going to be.

And in that case we get a curve that looks like this.

For voltages that are much lower than minus 85 millivolts, dV/dt is positive.

For voltages between minus 85 millivolts and around

plus 58 millivolts dV/dt is, is negative.

And then for voltages above minus 58 millivolts dV/dt is positive again.

So what happens when this crosses zero, right.

There are two, there are two voltages here

for which dV/dt instantaneously will be equal to zero.

And if dV/dt is equal to zero, that means that if you change voltage and you

put it at exactly that voltage right there

where it's not changing, then it's going to stay there.

And because these are voltages where the derivative of the voltage is zero.

These are, are what are known as fixed points.

So the fixed points in this case you can represent in this one dimensional example

by plotting with derivative on the y axis, tying the variable on the x axis.

And then seeing where the curve crosses zeri, right and so where

they, where, where this is zeri is the same where derivative is zero.

That means that if the voltage is at this level, and the derivative

at that, in that case is is zero, then it's going to stay there.

But then what we want to do is we want to say, well, what's going to happen

to the voltage when we move it away from one of these fixed points.

And this is where it gets more complicated and, and, and more interesting.

Right, if we start at minus 85 and we go negative.

This is a dV/dt is positive, that means

that voltage is going to move to the right.

Voltage is going to go up because the derivative is positive.

So when we change to minus 85 millivolts we have a positive dV/dt.

If we were to change something like minus

70, in this case the derivative is negative.

And so therefore we can draw the arrow this way.

And then if we were to go from minus 85 all the way to like minus

55, we would have a positive dV/dt, and so we would draw the arrow this way.

So that's why it's actually useful to, to draw this kind of plot, where you

have the derivative on the the y axis, and the variable on the x axis.

Is because knowing if the derivative is positive or negative in particular

ranges will tell you, well which way is the system going to evolve.

Which way are things going to change when we when we

move the variable, in this case, voltage to that particular location.

So here we have an arrow going to the

right, indicating that the volt, the derivative is positive.

Here we have an arrow moving to

the left, indicating that the derivative is negative.

And here we have an arrow moving to

the right again, indicating that the derivative is positive.

Why is it useful to plot the the arrows this way to show how the system

evolves when we change the when we change the variable, in this case the voltage.

Well this is a way that we can tell whether

our fix points are stable fix points or unstable fix points.

What do we mean by that?

Well, this fix point where the derivative

crosses zero has an arrow pointing towards it.

And an arrow pointing on the left, and an arrow pointing towards it on the right.

That means that we start at minus 85 and we move away from it, it's

going to go back to minus 85, because that's the way that it's pointed, right?

Similarly, if we go from minus 85 to minus 75, the arrow is

pointing negative, it's pointing to the left, it's pointing back toward minus 85.

And so, the voltage as we go back is going to go back to minus 85.

So, we can now delineate this fixed point as a stable fixed point because

deviations away from that fixed point will always push the system back towards it.

This fixed point at minus 58, in contrast, is an unstable fixed point.

What happens if you move negative to minus 58, well the derivative in that case

is going to be negative, and that means

that it's voltage is going to continue moving negative.

Similarly, if you're at minus 58 and you

get nudged a little bit positive, the derivative

is going to be positive and that means

it's going to continue to move away from that.

So, this fixed point has arrows pointing away from it.

That's how we identify it as an unstable fixed point.

And this fixed point has arrows pointing towards it.

That's how we identify it as a stable fixed point.

And we can perform this experiment numerically by changing voltage and then

allowing the system to evolve, then

integrating the equations numerically using Euler's Method.

And we see exactly what we can see over here on the left in schematic format.

We're starting at minus 85 millivolts.

When we hyper-polarize the cell, when we move the voltage down

to minus 95 millivolts, that's the black curve, what do we see?

It comes back to minus 85.

When we go up to minus 75, which is a Wrenn

curve, we see, we instantaneously change it and then it goes back.

Same thing with minus 65 which is the green curve.

it goes up to that instantaneously and it goes back.

And we can identify mi, minus 75 and minus 65.

On this, on this plot, and see that yes, the system

will go back to minus 85, when we're in that regime.

But if we start at minus 85, and we move the

cell all the way to minus 55, what do we see?

The voltage continues to go up.

And then it reaches a peak, and then it

continues to evolve, and that's because there are multiple

differential equations in this model.

There's somewhere around the order of

eight differential equations in this model.

But we can see that small deviations from minus

85, the system will come back to minus 85 millivolts.

But a large deviation, if you cross this unstable fixed point, then you're going to

move away from that unstable fixed point and move away from minus 85 millivolts.

And so by doing this numerical equation here, we can

confirm what we see graphically over here on the left.

That between mi, between minus 85 and minus 58, or something that's negative

to minus 85 is going to move back to the resting state, to minus 85.

So therefore that's a stable fixed point.

Small deviations away from minus 58 are going to cause action potentials or

return to the states and therefore this one is an unstable fixed point.

What we're going to do next is we're going to learn to

analyze these stable fixed points and unstable fixed points more rigorously.

The example I just showed of of a cardiac

action potential, we were plotting everything as the derivative

of voltage, minus voltage but it is in fact

a model that has eight differential equations in that case.

And, what's happening with eight differential equations

you can't, you cannot visualize everything graphically.

So now what we are going to do is we're going to

look at a two variable model where you can analyze things graphically.

And that's going to help us develop the tools we need

of dynamical systems, to be able to analyze these models.

And the two dimensional example we're going to use

describes the biological phenomenon of glycolytic oscillations in yeast.

And here's an, an example from a, from a pretty famous paper showing what happens.

There's an experimental procedure that you have to undertake in order

to get to, in order to observe these oscillations, where you condition

the cells, you starve them, and then when you add glucose, if

you look at oxygen consumption in this case, what do you see?

You see it goes up and down, up and down, up and down.

And these oscillations will continue for a very long time.

And this is a phenomenon that's been, that's

been known to occur for quite a while.

And in fact, many mathematical models of this process

have been developed and have been published over the years.

And, we're going to analyze one that was published

in, in the year 2000 by Bier et al,

and the citation for this for this model

that we're going to consider is given right here.

And this is the, the example system that we're going to use to introduce

concepts of stability and to show you how to analyze stability of ODE systems.

Now we want to ask what's the structure of the Bier

et al model and it's summarized in this diagram here.

The Bier et al model has two variables so the two variables that are computed that

are evolved with respect to time are glucose

inside the cell and ATP inside the cell.

Now we'll go through each step in the process.

Model simulates transport of glucose from the outside of the cell

to the inside of the cell through this rate called, called VN.

Then once glucose is inside the cell it can

get converted into ATP and this step here represents glycolysis.

in, in real life glycolysis repre, is occurs through the action

of several different enzymes, there are many different steps in this process.

But to keep things simple Bier et al represented lumped all of

these steps into, into one, and used the single rate constant k1.

The most important enzyme in this process, the rate limiting step, is

an enzyme that some of you may have heard of called phosphofructokinase.

So, this this rate k1, again it represents the

action of several steps in the process, several different enzymes.

But it can be generally considered primarily

the action of phosphofructokinase activity within the cells.

Once ATP is produced, it can get consumed.

And again, this is another place where, where Bier et al made an approximation.

There are several different ATPases that act,

that act within cells, each one of which

has, has different activity to keep things simple

Bier et al lumped all the ATPases together.

And then there's one more aspect of this model that's really

important and that's this dashed arrow, right, I drew down here.

Glycolysis is an interesting reaction in that it produces ATP

from glucose, but it also requires ATP to be initiated.

So when ATP goes up there are a couple of actions.

One is that ATP goes up, then ATP gets consumed, but the other thing

that happens is that this reaction here,

this production of ATP occurs more quickly.

And so we draw this as a regulatory arrow as ATP goes up.

Then you get more conversion of glucose into ATP.

So, we can look at the equations for the differential

equations for ATP and the differential equations for glucose and

we can understand the four terms here as the action

as either the steps that increase glucose or decrease glucose.

The steps that increase ATP or decrease ATP.

We'll look at glucose, dG dt first this ODE for glucose, initially.

Right, if you have more transport of glucose

into the cell, then glucose concentration's going to go up.

And then this term here represents conversion

of glucose into AT, ATP through glycolysis,

through the action of several different enzymes,

but phosphofructokinese being the most important one.

Now if we look at the things that

can either produce or consume ATP, this positive

term here represents conversion of glucose into ATP,

and again it depends on both glucose and ATP.

And then this term here represents consumption

of ATP through the action of ATPases.

And there are four four primary parameters in the Bier et

al model and the default values of three of them VN, K one and KP are given here.

And there is a fourth model, a

fourth model parameter that we haven't really discussed

yet, and that is the the KM for the action of the, of the ATPases.

And let's look at what happens in the, in the Bier et al model when we

change this fourth parameter, the Km for the

ATPases and the kls constants for the ATPases.

When Km is equal to 13, we see

these sustained oscillations of both glucose and ATP.

Glucose here is plotted in black.

ATP is plotted in red.

And, and you see that glucose goes up and goes

down, goes up and goes down, and this continues indefinitely.

Similarly, ATP exhibits oscillations, and they continue indefinitely.

What happens if we set Km equal to 20 instead of 13.

Well at the very beginning of this simulation, we

see oscillations, but what do we notice about these oscillations?

As time goes on the oscillations keep getting smaller, and smaller, and smaller.

And this is what we would refer to as,

as damped oscillations and with Km equal to 20.

And then finally what we notice is the

oscillations get so small you can't see them anymore.

So eventually glucose settles in at a constant

value and ATP settles in at a constant value.

So just like what we saw with the three component repressive network,

we can see qualitatively different behavior

as we change one of the parameters.

We see sustained oscillations over here, with Km equals 13, and we see

damped oscillations, and eventual settling at steady

levels of glucose and ATP over here.

So what we want to do next is we want to ask, how

can we understand the qualitatively different behavior

that we see in these two cases.

[SOUND] What we want to introduce now are

what we call phase plane techniques for analyzing two dimensional systems of ODEs.

And what we mean by phase plane is that

we're not going to plot glucose and ATP versus time.

Instead of plotting glucose versus time, and ATP

versus time, we want to plot glucose versus ATP.

So, glucose will be on the y axis in this case, and ATP will be on the x axis.

So let's take the time courses that we simulated on the last slide, for

Km equals 13 on the left, and for Km equals 20 on the right.

And plot them in the phase plane where ATP is

on the x axis and glucose is on the y axis.

And what we see with ATP what we see in

the phase plane with Km equals 13 is this loop.

And this loop makes sense intuitively because we saw that glucose

would oscillate, and ATP would oscillate and those would continue indefinitely.

And when you have two variables that are, that are both oscillating, what you see

when you plot one versus the other is

something like this that looks like a loop.

What about for case of Km equals 20

where you saw these damped oscillations, where you

saw these small fluctuations in the beginning but

eventually they settle down to a steady level.

Well, what this trajectory looks like in the phase plan is a something that

looks like this, it starts here and then it goes like this and it spirals.

And it goes down to some level and it eventually, it keeps spiraling and then

eventually it gets to some steady level where

ATP stays constant and glucose stays constant forever.

So the terminology we introduced to describe this is on the left.

We say that glucose and ATP oscillate indef, indefinitely

in something that we call a stable limit cycle.

In contrast on the right, we see glucose and ATP that converge to a stable fixed

point and in the next lecture we will

discuss these, this terminology in, in more detail.

Why do we think it's useful to plot these in the phase plane?

When you get similar information when you

plot glucose versus time and ATP versus time.

Well, one of the reasons why the the plotting things in the phase plane

is useful is that what you care about is which direction the system is moving.

And the direction is determined by the derivative of ATP with

respect to time and the derivative of glucose with respect to time.

And so at any given location in your phase

plane and any combination of ATP and glucose, you

can calculate the root of ATP with respect to

time and the derivative of glucose with respect to time.

And this defines a vector in the phase plane,

and this vector tells you how the system is

moving, how ATP is changing with respect to time

and how glucose is changing with respect to time.

So we can understand this directionality of our of

our two, of our two dimensional derivative vector, like this.

This is our differential equation for glucose,

this is our differential equation for ATP.

And what we saw was with Km equals 13 we, we saw the stable limit cycle.

One question becomes well, which way is it moving around this stable limit cycle?

Is it going clockwise?

Like I've drawn it here with the arrows,

or is it, would it be going counterclockwise?

Well, we can understand that by looking at these differential equations here.

What happens if we had very large values of glucose and, and ATP?

[SOUND] Glucose, the changing glucose with respect to time would be negative.

Because Vn in this case is a constant.

And glucose and ATP are multiplied together, in this case,

on this term, that has a negative in front of it.

So eventually, this term would become greater than this term Vn.

And glucose would be, would would be decreasing with respect to time.

We'd have a negative derivative.

So that's how we know that this arrow for glucose when we are up

here at the top right corner of this phase plane, should be pointed down.

What about ATP?

Well, when glucose and ATP both get large this term here, this first

term keeps getting bigger and bigger, this positive term gets bigger and bigger.

This term here, which has a negative in

front of it, doesn't continue to get bigger.

As ATP grows this term is eventually going to saturate.

And so this term is going to go up for low levels of ATP but

then it'll get to a point where it cannot continue to go up any more.

So clearly this term is going

to eventually overcome this negative term here.

So the change in ATP with respect to time for

large values of glucose and ATP is going to be positive.

So, just by looking at these differential equations qualitatively

we can draw this arrow up here for large levels

of glucose and ATP, and we know that it's pointing

to the right because the derivative of ATP is positive.

And we know that it's pointing down because

the derivative of glucose in this case is negative.

And so now that we know that this arrow is pointing this way,

we can then deduce how things are going through the rest of the system.

And we could have made the same sort of argument if

we had picked very low levels of glucose and ATP, for instance.

In the next lecture, we'll we'll talk about how we

can deduce when these arrows are going to switch directions.

Why does it go from pointing down and to the right

over here, to pointing down and to the left over here?

So, we'll talk about that in the next lecture.

But for now, we just want to introduce this concept of plotting,

of looking at systems and how they involve in the phase plane.

Which means you have one variable on the y

axis, and another variable here on the x axis.

To summarize this lecture, a fixed point of a dynamical system

is a set of variables where all derivatives are equal to zero.

A fixed point can be stable.

Which means that small perturbations away from that fixed point will

cause the system to evolve, and return back to that fixed point.

Conversely, as we saw in the one dimensional

example, we can also have unstable fixed points.

Which means that if you have a small perturbation away from an

unstable fixed point, it will continue to move away from that fixed point.

Finally, with two variable systems it's helpful

to plot one variable versus the other

and that's what we refer to as analyzing the system in the phase plane.

Lecture 8 - Introduction to Dynamical Systems - Part 3

系统生物学的动力学建模方法

与讲师见面