Lecture 30 - Stochastic Modeling - Part 2

Loading...

En provenance du cours de Icahn School of Medicine at Mount Sinai

Dynamical Modeling Methods for Systems Biology

123 notes

Icahn School of Medicine at Mount Sinai

Dynamical Modeling Methods for Systems Biology

123 notes

Cours 4 sur 6 dans Specialization Systems Biology and Biotechnology

An introduction to dynamical modeling techniques used in contemporary Systems Biology research. We take a case-based approach to teach contemporary mathematical modeling techniques. The course is appropriate for advanced undergraduates and beginning graduate students. Lectures provide biological background and describe the development of both classical mathematical models and more recent representations of biological processes. The course will be useful for students who plan to use experimental techniques as their approach in the laboratory and employ computational modeling as a tool to draw deeper understanding of experiments. The course should also be valuable as an introductory overview for students planning to conduct original research in modeling biological systems. This course focuses on dynamical modeling techniques used in Systems Biology research. These techniques are based on biological mechanisms, and simulations with these models generate predictions that can subsequently be tested experimentally. These testable predictions frequently provide novel insight into biological processes. The approaches taught here can be grouped into the following categories: 1) ordinary differential equation-based models, 2) partial differential equation-based models, and 3) stochastic models.

À partir de la leçon

Stochastic Modeling

The description goes here

Lecture 29 - Stochastic Modeling - Part 1 16:22

Lecture 30 - Stochastic Modeling - Part 2 19:29

Rencontrer les enseignants

Eric Sobie, PhD
Associate Professor
Department of Pharmacology and Systems Therapeutics, Systems Biology Center New York (SBCNY)

Hello.

This is a second lecture on methods for solving Stochastic Mathematical Models.

In the second part of the lecture, we're going to discuss practical aspects.

How can we solve stochastic mathematical models?

In the first part,

we talked about what is, what we meant by a stochastic mathematical model.

And we talked about how conceptually

one has to think of things in terms of probabilities of reactions occurring,

rather than rates at which reactions occur.

Now we want to look at it from a practical point of view.

How do we solve stochastic mathematical models?

I'm going to introduce a couple of different approaches for

solving stochastic mathematical models.

First a simple brute, what you might call a brute force method for

solving these stochastic models.

I call this a stochastic Euler's method, because it's analogous to what we saw with

Euler's method for solving ordinary differential equations.

And then we want to introduce a more sophisticated one

which is known as Gillespie's algorithm.

And then we're going to finish this,

this lecture with that with a description of Gillespie's algorithm.

And we're going to show some MATLAB code in order to solve these stochastic models.

And the application that we're going to use is stochasticating of

the Hodgkin-Huxley potassium channel.

In the previous lectures on the Hodgkin-Huxley model of

the neuronal action potential, we discussed in detail

how those investigators developed their their models for the potassium channel.

And here were going to show how we, one can implement a stochastic version of

that, that incorporates the random opening and closing of individual channels.

Well start be reviewing one of the most important concepts that we discussed in

the last lecture.

Remember we talked about the fact that the,

the key in stochastic models is to convert from rates that reactions occur

into probabilities that reactions might occur.

Where a rate which we are used to thinking of with ordinary differential equation

models is expressed as the number of transitions per unit time, for

instance, reactions per second.

Now for stochastic models, we need to think in terms of probabilities.

Meaning, what's the chance that a reaction is going to occur

in a given amount of time?

And in the last lecture we went through this procedure by which one can convert

from rates into probabilities.

And we're not going to go over this in details again, but we had two steps here,

one is to convert from concentration per time to molecules per time.

And then the second step would be to convert from molecules per time

into a probability that an event will occur in a given time interval, and

this procedure here is going to be essential

in our development of stochastic mathematical models.

The other important concept I want to review goes back to our lectures from

the Hodgkin-Huxley model of, of the action potential, the ordinary differential

equation model, which is how Hodgkin and Huxley developed their model of

potassium conductance, which then lead to their model of potassium current.

We had these traces here which were obtained under voltage clamp

conditions for pota, potassium conductance as

a function of time at many different voltage clamp steps.

And one of the thing we, one of the things we discussed when we talked about

the Hodgkin-Huxley model is that those investigators noticed that changing

the voltage changes both the steady-state level of the potassium conductance, and

the rate at which it goes up.

In other words, it goes up faster at plus 60 millivolts than it does at these other

these other voltages.

They also noticed that the time course of potassium conductance increased

is similar to an exponential raised to a power, and

they specifically thought it looked like the fourth power.

So this is what happens, this blue curve is what you have if you have a,

exponential raised to the first power,

the red curve is an exponential raised to the second power,

and then the magenta curve here is this exponential raised to the fourth power.

And Hodgkin and Huxley deduced that these facts suggested the following model where

they introduced a variable n, which is a variable between 0 and 1,

representing the fraction of particles in our, that are in a permissive state.

And they said that the conductance of the,

of potassium conductance was proportional to n raised to the fourth power.

In other words you, if you had you had four particles and

all four of them had to be in the permissive state in order for

the the membrane to be allowing potassium current to go through.

So, in their model potassium conductance was equal to some maximum

potassium conductance times n raised to the fourth power.

And they simulated the transitions of between the non-permissive

state 1 minus n and the permissive state n, with these variables alpha and

beta, where in general alpha and beta are functions of voltage.

And then remember what I, I said just a minute ago, is that your, your variable n,

the gating variable, because it represents a fraction of particles in the permissive

state, this will always be some number between 0 and 1.

This is just to review our what we learned before about potassium con,

conductance model from Hodgkin and Huxley, so that we can introduce a stochastic

model for the gating of a single potassium channel.

Because in the previous slide we had the potassium conductance was proportional to

n raised to the fourth power,

this implies the following model of a single potassium channel.

It looks something like this, which looks pretty complicated, but

as we'll see it's not nearly as complicated as you might think.

So you can't just say that the channel can be closed or open.

We have to posit several different closed states here which we call C0,

closed 0, closed 1, closed 2, closed 3, and then open state.

Why do we have to introduce all these states here?

Well, the idea here, and the reason we number them as we did 0, 1,

2, and 3, is that this state over here, C0,

is when we have zero out of four of our particles in the permissive state.

This is for one out of four permissive.

This is for two out of four permissive, et cetera, and then as we just discussed in

the previous slide, when all four particles are in the permissive state then

that is when the channels open, and that is a consequence of this equation here,

where the conductance is proportional to n raised to the fourth power.

In other words, you have four, gates on your sodium channel.

All four of them have to move into the permissive state in order for

the channel to open.

Now, what are the rates that govern the transition from C0 to C1,

from C1 to C2, et cetera?

It's actually something that we can get just based

on knowing how each particle transitions from permissive to non-permissive.

This is what we had in the deterministic Hodgkin-Huxley model, where the transition

from the non-permissive state, 1 minus n, to the permissive state n, occurs at a,

a rate alpha, and the transition in the other direction occurs at a rate beta.

So, what that tells us is that, when we're transitioning from C0 to C1,

the overall rate is 4 times alpha.

The reason we have a 4 in front of it is because when zero out of four are perm,

are permissive, that means that there,

there are four particles that might be able to transition.

That's why we multiply this one by 4.

Here there are three particles that might be able to transition.

So it's 3 alpha, and then 2 alpha, and

then 1 alpha, and then it's the opposite going back in the other direction.

When all four particles are in the permissive state, then you have four

particles that could possibly transition back into the non-permissive state.

So that's why we multiply beta by 4, and then over here we multiply beta by 1.

Now that we've established that this is our model of the channel, converting this

into a stochastic gating model in this case is relatively straightforward.

If our model of a single potassium channel looks like this,

as we just discussed in the last slide, we can write down a simple procedure for

stochastically simulating channel gating that looks like this.

At any given time, the channel's going to be in a particular state.

Let's just assume for the sake of argument that,

at this moment, the channel is in state C1.

Extending this to all the other states is, is straightforward but, at any given time,

you're going to be in a particular state.

Let's say that we're in C1.

What you do in your algorithm next is you compute the probabilities that you, well,

are going to transition out of that state.

Now the probability that you're going to go

from C1 to C2 is equal to 3 times alpha times your times step.

Remember, this is what we discussed in the previous lecture,

that the probability you're going to have a transition in a particular amount of

time is whatever that that propensity is times your time step.

So, the probability you're going to go from C1 to C2 is dt times 3 alpha.

The probability that you're going to go from C1 to C0 is equal to dt times beta.

Then you pick some uniformly distributed random numbers.

If a uniformly distributed random number is less than the transition probability,

then the channel is going to change states.

In other words, if you pick a random number that's less than dt times 3 alpha,

that's going to give you a transition from C1 to C2.

If you pick a uniformly distributed random number that's less than dt times beta,

you're going to transition from C1 to C0.

Now you might be saying to yourself, well, what if you pick

a random number that will tr, tell it's going to transition both to C2 and to C0?

That's that's obviously impossible for both of those to occur in dt,

and that tells you something about what sort of time step you have to choose for

an algorithm like this.

You need to choose dt very, very small, so that the probability of

both transitions occurring has to be, become vanishingly small.

And so I think the rule of thumb that a lot of people use for

these sorts of these algorithms is you want the probability of any

single transition to occur, occurring to be much less than 1%.

Therefore, the probability that multiple transitions will be predicted to occur

is very, very, very small.

Once you do this, you update the state if necessary.

So, you pick a random number.

You say, okay maybe, I'm going to have a tra, transition from C1 to C2 or maybe,

I'm going to have a transition from C1 to C0.

Most of the time I'm not going to have either transition.

You update it if you need to, you record your results, and

then you move on to the next time stop.

This is actually a rather straight forward algorithm for

simulating the stochastic gating of this potassium channel.

And we can see what results we get when we implement this

relatively simple stochastic model of potassium channels.

We're going to simulate a voltage clamp step from minus 60 millivolts

to 0 millivolts.

Remember, we discussed these sorts of experiments in our

lectures on the Hodgkin-Huxley model where you hold the axonal membrane

at minus 60 millivolts, and then you depolarize to 0 millivolts, and

the potassium channels will open.

In this case it's not enough just to run the simulation once.

You have to run it multiple times over and over again, because when you run

the simulation on once channel you'll see the single channel opening and

closing randomly.

And this is what the output looks like when you run multiple trials with

a single channel.

I think the depolarization occurs right around here maybe after one millisecond or

so, and sometimes there's a long delay before it opens,

other times there's a short delay.

You can also see that how long it's open for is randomly distributed.

Sometimes it's open for a long time, sometimes it's only open for

a short time and occasionally you might not get any openings at all in,

in this time interval, that's what we see with this cyan one here.

So what each one of these traces is, is showing is either the closed state,

which is the lower level, or the open state which is the upper level,

and by convention we denote this as little ik it was,

it was a big I for, for current though all the potassium channels.

And little ik tells us, tells us the current through one channel, and

we can see individual channels in this case are opening and closing at random.

We can also compare the results of this model with what we get

from the deterministic or ordinary differential equation model.

If we average 100 of these voltage clamp steps,

we get something that looks like this, an increase in

potassium current is a function of time with a lot of noise on it.

And the noise reflects the, the randomness that the channels open and

close at random, and but

this looks like very much like what we see with the ordinary differential equation

model solution, where the the current goes up smoothly as a function of time.

So as our number of of trials with a single channel gets bigger,

and bigger, and bigger the solution with

our stochastic solution is going to approach our deterministic solution and

it has to because we're essentially simulating the same model here.

It's just in one case we're, we're incorporating some randomness,

in the other case we're not.

In the last slide, we discussed a very simple algorithm for solving stochastic

mathematical models, what I like to call the stochastic Euler's method.

But one important aspect of this very simple algorithm

is that during most intervals, no transitions occur.

In the examples I showed on the last slide of the the, potassium channel that was

opening and closing stochastically, the time step I used was four microseconds.

And during most of these four microsecond intervals if

the channel starts in a particular state at the end of that four microsecond

interval it's going to remain in that state.

So then we can ask, well what if, instead of calculating for

each time step whether or not a random event occurs,

what if instead we computed the interval until the next random event?

Now the thing about how might implement an algorithm like

this we need to remember some things from probability theory.

If stochastic transitions occur on average n times per second then

the average interval between transitions is going to be 1 over n seconds.

In other words, if our transitions are occurring on average ten times per second,

then our average interval is going to be one-tenths of a second or 100 milliseconds

between events, and that's required to get ten transitions on average per second.

Furthermore, we also know that these intervals are going to be exponentially

distributed, and what I mean by that is if you collected lots of intervals,

one after the other, after the other, after the other and

you plotted a histogram of these intervals, it would look like this,

where very, very long intervals would become more and more rare.

Short intervals would be more common and the mean value you would have

would be equal to 1 over n, and it would be indicated something like this.

And the histogram going from like, the shorter intervals to the longer intervals

would follow the, the shape of an exponential.

And so, the intervals in this case would be exponentially distributed

with a mean value of 1 over n in this case.

And this fundamental principle of probability theory forms the basis for

the algorithm we're going to discuss next which is called Gillespie's algorithm.

Why do we call this Gillespie's algorithm?

Well, as you can probably guess,

it's because it was invented by a man named Gillespie.

And this is probably still the most commonly used stochastic algorithm for

solving these sorts of systems.

It comes from this paper here, published all the way back in 1977.

That's a a mark that, that you've published something that's very

influential, that's an important piece of work when people actually,

still talking about it all these years down the, down the line.

So, Daniel T.Gillespie published this paper in 1977,

in the Journal of Physical Chemistry, that described this,

this algorithm that people still commonly use these days.

Gillespie's algorithm is similar to the stochastic

Euler's method that we discussed before, but it's, it's a little bit more complex,

but in return for this complexity it's a lot faster.

Gillespie's algorithm involves selecting two random variables at each time step.

And I'll go through the process, so we can see how how this works, and

how this algorithm can be quite efficient.

First thing you do is you compute the overall propensity for

any reaction to occur.

We're going to go back to that same example with the the stochastic gating of

the potassium channel, and we're going to assume that

at the current time step the potassium channel is in state C1.

And if it's in state C1, there are two transitions that are possible, right?

It can either go from C1 to C2, which has a propensity of 3 alpha, or

it can go from C1 back to C0, which has a propensity of beta.

So the overall propensity is the sum of these two, 3 alpha plus beta.

The next thing you do with Gillespie's algorithm is you select a random number

from an exponential distribu, distribution with mean of 1 over the, the propensity.

And what this random number from the exponential distribution tells you

is this is the interval until the next reaction occurs.

And then the third thing you have to do is select ano, another random number.

In this case it's a uniformly distributed random number which we're

going to call little r, and little r will tell you which reaction occurs.

So in this case, you start in state C1, you can either move to C2 or you can

move to C0 and we discussed the overall propensity for a reaction to occur.

We select the the interval that tells us when the next reaction is going to occur,

but we don't know if it's going to be moving from left to right, or

moving from right to left in this case.

The second uniformly-distributed random number, little r,

tells us which reaction occurs.

And in this case, it depends on the relative

magnitudes of this propensity beta and this propensity 3 alpha.

So if little r is less than beta over beta plus 3 alpha,

that tells us we have a C1 to C0 transition.

On the other hand if it's between this ratio, beta over

beta plus 3 alpha and 1, that tells us we have to have a C1 to C2 transition.

Then the fourth thing we do, is we update the system state,

repeat move onto the next iteration of Gillespie's Algorithm.

In summary then, what we've seen is that stochastic mathematical

models are simulated with algorithms that generate random numbers, and these random

numbers tell you whether a reaction occurs or when the reaction occurs.

Two examples of algorithms are used for stochastic systems are the stochastic

Euler's method, and this one involves a fixed time step, and

Gillespie's algorithm which involves a variable time step.

Finally, we've also seen that a simple stochastic model of the neuronal potassium

channel can reproduce the average behavior that's seen in the Hodgkin-Huxley

potassium current model.

Explorer notre catalogue

Rejoignez-nous gratuitement et obtenez des recommendations, des mises à jour et des offres personnalisées.

Coursera propose un accès universel à la meilleure formation au monde, en partenariat avec des universités et des organisations du plus haut niveau, pour proposer des cours en ligne.

© 2019 Coursera Inc. Tous droits réservés.

Télécharger dans l'App Store

Disponible sur Google Play