The following content is provided under a Creative
Commons license.
Your support will help MIT OpenCourseWare
continue to offer high quality educational resources for free.
To make a donation or view additional materials
from hundreds of MIT courses, visit MIT OpenCourseWare
at ocw.mit.edu.
PROFESSOR: Today what we want to do
is discuss various approaches that you
might want to take towards trying to understand
stochastic systems.
In particular, how is it that we might model or simulate
a stochastic system?
Now, we will kind of continue our discussion of the master
equation from last time.
Hopefully now you've kind of thought about it a bit more
in the context of the reading.
And we'll discuss kind of what it means to be using the master
equation and how to formulate the master equation for more
complicated situations, for example, when you have
more than one chemical species.
And then we'll talk about the idea of this Gillespie method,
which is an exact way to simulate stochastic systems,
and it's both exact and computationally tractable as
compared to what you might call various naive methods.
And the Gillespie method is really sort of
qualitatively different from the master equation
because in the master equation, you're
looking at the evolution of probability distributions
across the system, whereas the Gillespie method is really
a way to generate individual stochastic trajectories.
So if you start with somehow similar initial conditions,
then you can actually get-- you can
get, for example, the probability distributions
from the Gillespie method by running
many individual trajectories.
But it's kind of conceptually rather different
because of this notion of whether you think
about probabilities or you're thinking
about individual instantiations of some stochastic trajectory.
So we'll try to make sense of when you might
want to use one or the other.
And then finally we'll talk about this Fokker-Planck
approximation, which, as the reading indicated,
for intermediate ends, it's useful to make
this kind of continuous approximation,
and then you can get a lot of intuition
from your knowledge about diffusion
on effective [INAUDIBLE] landscapes.
Are there any questions about this or administrative things
before we get going?
I just want to remind you that the midterm is indeed
next Thursday evening, 7-9 PM.
If you have a problem with that time,
then you should have emailed [? Sarab. ?]
And if you haven't emailed him yet,
you should do it right now.
And-- yes.
All right.
So let's think about the master equation a little bit more.
Now before what we did is we thought about the simplest
possible case of the master equation, which
is, if you just have something being created
at a constant rate and then being degraded at a rate that's
proportional to the number of that chemical species.
And I'm going to be using the nomenclature that's
a little bit closer to what was in your reading, just
for, hopefully, clarity.
And I think that some of my choices from last lecture
were maybe unfortunate.
So here, this is, for example, m would be the number of mRNA,
for example, in the cell.
This is the rate of creation of the mRNA,
and then the rate of degradation of the mRNA.
So m is the number of mRNA.
And if we want understand gene expression,
we might include an equation for the protein,
so we might have some p dot, where some Kp.
Now-- oh, sorry.
Again, I always do this.
All right.
So we're going to have this be an n dot.
So now n is going to be the number of the protein.
Now this really is kind of the simplest possible model
that you might write down for gene expression that includes
the mRNA and the protein.
So there's no autoregulation of any sort.
It's just that the mRNA is involved
in increasing the protein, but then we
have degradation of the protein as well.
So what we want to do is kind of try
to understand how to formulate the master equation here.
But then also, we want to make sure
that we understand what the master equation is actually
telling us and how it might be used.
So first of all, in this model, I want to know
is there, in principle, protein bursts?
So before we talked about the fact that in-- at least
in [? Sunny's ?] paper that we read--
they could observe protein bursts,
at least in those experiments in e Coli.
Question is, should this model somehow exhibit protein bursts,
and why or why not?
I just want to see where we are on this.
I think this is something that, depending
on how you interpret the question,
you might decide the answer is yes or no.
But I'm curious-- I think it's worth discussing
what the implications are here.
And the relevant part of this is going
to be the discussion afterwards, so I'd
say don't worry too much about what you think right now.
But I'm just curious.
This model, does it include, somehow, protein bursts?
Ready?
Three, two, one.
OK.
So we got-- I'd say at least a majority of people
are saying no.
But then some people are saying yes.
So can somebody volunteer why or why not?
Yes?
AUDIENCE: I think the difference is if we're-- are we using this
in a continuous fashion or are we using this in a discrete
fashion [INAUDIBLE].
PROFESSOR: Yeah.
OK.
All right.
So he's answered both possible sides of the argument.
And the point here is that if you just
simulate this from the standpoint--
certainly, for example, this continuous,
this discrete-- so if you just simulate
this as a deterministic pair of differential equations,
then will there be bursts?
No.
Because everything is well-behaved here.
On the other hand, if we go and we do a full Gillespie
simulation of this pair of equations, then
in the proper parameter regime, we actually
will get protein bursts, which is, in some ways,
weird, that depending upon the framework that you're
going to be analyzing this in, you
can get qualitatively different behaviors for things.
But there's a sense here that the deterministic, continuous
evolution of these quantities would be the average over many
of these stochastic trajectories,
and the stochastic ones do have bursts,
but if you average over many, many of them,
then you end up getting some well-behaved pair of equations.
So we'll kind of try to make sense of this more later on.
But I think this just highlights that you
can get really qualitatively different behaviors
for the same set of equations depending
upon what you're looking at.
And these protein bursts can be dramatic events, right,
where the protein number pops up by a lot.
So this really, then, if you look
at the individual trajectories here,
they would look very different whether you
were doing kind of a stochastic treatment
or the deterministic one.
Can somebody remind us the situation
in which we get protein bursts in the stochastic model?
In particular, will we always get these discrete protein
bursts?
Or what determines the size of a protein burst?
Yes.
AUDIENCE: Does it have to do with the lag
time between when the mRNA is created [INAUDIBLE]?
PROFESSOR: OK.
Right.
So there's a lag time between the time that mRNA is created,
and then the next thing would be--
AUDIENCE: When the protein [INAUDIBLE].
PROFESSOR: When the protein is [? totaled-- ?] OK.
So there are multiple time scales, right?
So after an mRNA is created, and that's
through this process here-- so now out pops an mRNA--
now there are multiple time scales.
There's the time scale for mRNA degradation.
That goes as 1 over gamma m.
There's a time scale for protein degradation
after a protein is made.
That goes as 1 over gamma p.
But then there's also a time scale associated
with kind of the rate of protein production
from each of those mRNAs, and that's determined by Kp.
So we get big protein bursts if what?
What determines the size of these protein bursts?
Yes.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Right.
It's always confusing.
We talk about times.
But in particular, we have protein bursts
in the stochastic situation if we do a stochastic simulation.
And that's in the regime if Kp, the rate
of protein synthesis from the mRNA
is somehow much larger than this gamma m.
Have I screwed up?
Yes.
AUDIENCE: So this is also-- in the sense
of being different from the deterministic equations,
we probably also want the total number of mRNAs [INAUDIBLE].
Is that sort of--
PROFESSOR: Right.
Yeah, I think that it-- and the question of what
mRNA number you need.
I mean, it depends on what you mean by protein bursts.
I would say that so long as this is true,
what that means is that each mRNA will, indeed, kind of lead
to a burst of proteins being made, where the burst is,
again, geometrically distributed with some--
now there's another question, which is, are those protein
bursts kind of large compared to the steady state protein
concentration?
And that's going to depend upon Km and gamma p as well.
Is that--
AUDIENCE: Yeah.
So I guess [INAUDIBLE] which is, I guess it would also
depend on how big [INAUDIBLE].
PROFESSOR: All right, well-- and you're saying time resolution
in terms of just measuring--
AUDIENCE: Yeah. [INAUDIBLE]
PROFESSOR: Yeah.
Well, OK, but right now we're kind
of imagining that we live in this perfect world
where we know at every moment of time
exactly how many of everything there is.
So in some ways we haven't really said anything yet
about time resolution.
We're assuming that our time resolution and our number
resolution is actually perfect.
But still, depending upon the regime that you're in,
the protein numbers could look something
like-- so if you look at the protein number, which
is defined as this n as a function of time, then
in one regime, you're going to see where it's kind of low.
You get a big burst and then it kind of comes down, and then
a big burst, and then it kind of comes down, and burst,
and it kind of comes down, right?
So this is in the regime where you have infrequent mRNAs being
produced, and then large size bursts from each mRNA.
And then you kind of get this effective degradation
or dilution of the protein numbers over time.
And this distribution, if you take a histogram of it,
is what?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Right.
So I'm imagining that we look at this for a long period of time.
And then we come over here and we histogram it.
So now we come over here, we turn to the left,
we say number has a function of-- this is number n.
The frequency that we observe, some number of proteins.
So frequency.
And this is going to do something.
So what about-- it may not be a beautiful drawing,
but you're supposed to know the answer.
I'm trying to review things for you because I
hear that you have a big exam coming up,
and I want to make sure that--
Gamma.
It's a gamma, right?
So this is what we learned earlier.
So this is a gamma distribution.
And you should know what this gamma distribution looks like.
In particular, there are these two parameters
that describe this gamma distribution
as a function of underlying parameters in the model.
AUDIENCE: [INAUDIBLE]
PROFESSOR: Maybe.
I don't want to get too much into this because, well,
on Thursday we spent a long time talking about it.
Once we get going, we'll spend another long time
taling about it again.
But you should review your notes from Thursday before the exam.
So this thing is gamma distributed.
And if we looked at the mRNA number as a function of time
and we did a histogram of that, the mRNA distribution
would be what?
It's poisson.
So it's important to remember that just because I tell you
that a protein number is gamma distributed,
that doesn't immediately tell you exactly what you should
be expecting for the distribution of, say,
the number of protein as a function of time.
I mean, there are many different things
I could plot over here that would all kind of come down
to a gamma distribution over here.
So it's important to kind of keep
in mind the different representations that you might
want to think about the data.
So what we want to do now is we want to think a little bit more
about this master equation in the context of if we're
going to divide it up into these states.
Now I would say that any time that you
are asked to write down the master equation for something--
so now how many equations will the master equation-- I
say master equation, but there is really more than one, maybe.
So how many equations will be involved in the master equation
kind of description of this model?
Infinitely many.
But there were infinitely many already when
we had just one, when we just had the mRNA distribution.
Well, you know, infinite times infinite is still infinite.
So long as it's a countably infinite number.
But yeah, but it's still infinite, always.
All right.
So what we want to do is divide up the states.
So when somebody asks you for-- the equation's describing how
those probabilities are going to vary,
really what we're interested in is some derivative with respect
to time of some probabilities described by m,n.
We want to know the derivative with respect to time for all
m,n's.
So that's why there are infinite number, because m goes in one
direction, n goes in another.
Lots of them, OK?
Now it's always tempting to just write down this derivative
and then just write down the equation.
If you do that, that's fine, but I
would recommend that in general what
you do is you try to write a little chart out
to keep track of what directions things can go.
So for example, here we have the probability of being the m,n
state.
Now there's going to be ways of going here.
And this is going to be going probability of being an m plus
1,n.
What I'm going to do is I'm going to give you
just a couple minutes.
And in two minutes, I want you to try to write down as many
of the rates, the f's and n's that correspond
to all these transitions.
You may not be able to get through all of them,
but if you don't try to figure out some of them,
then you're going to have trouble doing it
at a later date.
Do you understand what I'm asking you to do?
So next to each one of these arrows,
you should write something.
So I'll give you two minutes to kind of do
your best of writing these things down.
All right.
Why don't we reconvene, and we'll see how we are?
So this is very similar to what we did on Thursday.
We have to remember that m's are the mRNAs,
and this is what we solved before,
where it's just a long row.
Now first of all, the mRNA distributions and the rates,
do they depend on the protein numbers?
No.
So what that mean about, say, this arrow as
compared to the arrow that would be down here?
It's going to be the same, because n does not
appear in that equation describing mRNA.
If we had autoregulation of some sort, then it would.
So let's go through.
All right.
What we're going to do is we're going
to do a verbal yelling out.
OK, ready.
This arrow.
AUDIENCE: Km.
PROFESSOR: This one here is, 3,2,1--
AUDIENCE: Km.
PROFESSOR: Km.
All right.
Ready, 3, 2, 1.
AUDIENCE: Gamma m times m.
PROFESSOR: Gamma m times m.
3, 2, 1.
AUDIENCE: Gamma n times m plus 1.
PROFESSOR: Gamma m times m plus 1.
Now remember that there are more mRNA over here then there
are here, which means that the rate of degradation
will increase.
Now coming here, now this is talking
about the creation and destruction
of the proteins, changes in n.
All right, this arrow here.
Ready, 3, 2, 1.
AUDIENCE: Kp times m.
PROFESSOR: It's Kp times m.
So this is the rate of creation, going from n minus 1 to n.
That's fine.
You know, I was looking at my notes from last year,
and I got one of these things incorrect, so-- and then,
OK, ready.
This one here, 3, 2, 1.
Kp times m.
So here the same rate, and should we be surprised by that?
So the number of proteins are changing,
but here it's the number of mRNA that
matters, because we're talking about the rate of translation,
right?
Now this one here, 3, 2, 1.
Gamma p times n.
And here, 3, 2, 1.
AUDIENCE: Gamma p times n plus 1.
PROFESSOR: Gamma p times n plus 1.
All right.
Perfect.
Now this is, of course, as you can imagine, the simplest
possible kind of set of equations that we
could have written down.
If you have other crazy things, you
get different distributions, if you have autoregulation
or if you have interactions of something
with something else, or the same thing, so forth.
But I think it's really very useful
to kind of write this thing down to clarify your thinking
in these problems.
And then you can fill out-- for change of probability,
you have mn.
You come here and you just go around
and you count, take all the arrows coming in,
and those are ways of increasing your probability.
And ways going out are ways of decreasing your probability.
Now in all those cases you have to multiply these raw rates
by the probabilities of being in all these other states.
So can you use the master equation
to get these probabilities if you're out of equilibrium,
out of steady state?
So that's a question.
So the master equation useful out of steady state?
Yes.
Ready.
3, 2, 1.
All right.
So we got a fair number of-- there is some disagreement,
but yeah.
So it actually-- the answer is yes.
And that's because you can start with any distribution
of probabilities across all the states that you'd like.
It could be that all of the probabilities at one state.
It could be however you like.
And the master equation tells you
about how that probability distribution
will change over time.
Now if you let that run forever, then you
come to some equilibrium steady state.
And that's a very interesting quantity,
is the steady state distribution of these probabilities.
But you can actually calculate from any initial distribution
of probabilities evolving to any later time t what
the probability would be later.
This comes to another question here.
All right.
So let's imagine that at time t equal to 0,
I tell you that there are m not mRNA and P
not-- I always do this.
I don't know, somehow my brain does not like this.
Because the P's we want to be probabilities.
We start with m not mRNA, n not protein.
And maybe it's a complicated situation.
We can't calculate this analytically.
So what we do is we go to our computer,
and we have it solve how this probability distribution will
evolve so that time T equal to some time--
if we'd like we can say this is T1.
I'll tell you, oh, the probability of having m and n
mRNA and protein is going to be equal to something P1.
Now the question is, let's say I then go
and I do this simulation again.
Now I calculate some other at time T1 again,
the probability that you're in the m,n state.
The question is, will you again get P1?
So this is a question mark.
And A is yes, B is no.
All right.
I'm going to give you 15 seconds.
I think this is very important that you understand
what the master equation is doing and what it is not doing.
AUDIENCE: [INAUDIBLE]
PROFESSOR: I'm sorry, what's that?
Right.
OK.
So I mean, this is just-- you know,
you program in your computer to use the master
equation to solve how the probabilities are
going to evolve.
I'm just telling you, start with some initial distribution.
And if you do it once, it says, oh,
the probability that you're going
to have m-- this time you're going to have mRNA proteins is
going to be P1, so it's 10%.
Great.
Now I'm asking just, if you go back and do it again,
will you again get 10%, or is this output stochastic?
It's OK that if you're confused by this distinction.
I think that it's easy to get confused by,
which is why I'm doing this.
But let's just see where we are.
Ready?
3, 2, 1.
All right.
So I'd say a majority again.
We're kind of at the 80-20, 75-25.
A majority here are saying that, yes, you
will get the same probability.
And this is very important that we understand
kind of where this where the stochasticity is somehow
embedded in these different representations
of these modelings.
The master equation is a set of differential equations telling
you about how the probabilities change over time
given some initial conditions.
Now we're using these things to calculate
the evolution of some random process,
but the probabilities themselves evolve deterministically.
So what that means is that although these things are
probabilities, if you start somewhere
and you use the master equation to solve,
you get the same thing every time you do it.
Now this is not true for the Gillespie simulation,
because that, you're looking at an individual trajectory.
An individual trajectory, then the stochasticity
is embedded in that trajectory itself,
whereas in the master equation, the stochasticity arises
because these are probabilities that are calculating,
so any individual instantiation will be probabilistic
because you are sampling from those different probability
distributions.
Now this is, I think, a sufficiently important point
that if there are questions about it,
we should talk about it.
Yeah.
AUDIENCE: How do you make the simulations?
Would you essentially-- can you take
a sum over different Gillespie?
PROFESSOR: So it's true that you can do
a sum over different Gillespie.
But we haven't yet told you about,
what the Gillespie algorithm is, so I can't use that.
But indeed, you can just use a standard solver
of differential equations.
So whatever program you use is going
to have some way of doing this.
And once you've written down these equations,
the fact that these are actually probabilities doesn't matter.
So those could have been something else.
So this could be the number of eggs, whatever, right?
So once you've gotten the equations,
then equations just tell you how the problems
are going to change over time.
Yeah.
AUDIENCE: Maybe this is a silly question,
but in practice, do you have to assume all the probabilities
are 0 above some number?
PROFESSOR: Oh no, that's not at all a silly question, because--
AUDIENCE: [INAUDIBLE]
PROFESSOR: Exactly.
Right.
And yes, it's a very good question.
So I told you this is an infinite set of differential
equations.
But at the same time I told you this master equation's supposed
to be useful for something, and kind of at the face of it,
these are incompatible ideas.
And the basic answer is that you have
to include all the states where there
is a sort of non-negligible probability.
We could be concrete, though.
So let's imagine that I tell you we want to look at the mRNA
number here.
And I tell you that OK, Km is equal to--
well, let me make sure.
Gamma m.
What are typical lifetimes of mRNAs in bacteria again?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Right.
Order a minute.
So that means that-- let's say this is 0.5 minutes minus 1.
To get a lifetime of around 2 minutes.
And then let's imagine that this is then 50 per minute.
So an mRNA is kind of made once a minute.
There's 50 of them.
That's a lot, but whatever.
There are a few genes.
Minute.
I wanted the number to be something.
So there's a fair rate of mRNA production.
Now how many equations do you think
you might need to simulate?
So we'll think about this.
First of all, does it depend upon the initial conditions
or not?
AUDIENCE: Maybe.
PROFESSOR: Yeah.
It does.
So be careful.
But let's say that I tell you that we start with 50 mRNA.
The question is, how many equations
do you think you might have to write down?
And let's say we want to understand this once it
gets to, say, the equilibrium.
All right.
Number of equations.
Give me a moment to come up with some reasonable options.
Well, these are-- let's say that this
could show up on your homework.
So the question is, how many equations
are you going to program into your intersimulation?
And it may be-- doesn't have to be exactly any of these number,
but order.
Do you guys understand the question?
So we need a different equation for each
of these probabilities.
So in principle we have-- the master equation
gives us an infinite number of equations.
So we have d the probability of having 0 mRNA with respect
to time.
That's going to be-- any idea what this is going to be?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Right.
So we have a minus Km times what?
Times p0, right.
So this is because if we start out down here at P0.
Now we have Km.
So I was just about to violate my rule
and just write down an equation without drawing this thing.
So it's Km times p0.
That's a way you lose probability,
but you can also gain probability
at a rate that goes as gamma m times P1.
So that's how this probability is going to change over time.
But we have a different equation for you for p1, for p2, for p3,
for p4, all the way, in principle, to p1,683,000,
bla bla bla, right?
So that's problematic, because if we
have to actually in our program code up 100,000,000 equations,
or it could be worse.
Then we're going have trouble with our computers.
So you always have to have some notion
of what you should be doing.
And this also highlights that it's really important
to have some intuitive notion of what's going on in your system
before you go and you start programming,
because in that case-- well, you're
likely to write down something that's wrong.
You won't know if you have the right answer,
and you might do something that doesn't make any sense.
So you have to have some notion of what
the system should look like before you even
start coding it.
My question is, how many of these equations
should we simulate?
OK.
Let's just see where we are.
Ready.
3, 2, 1.
OK.
So I'd say that we have, it's basically
between C and D. Yeah.
I would say some people are maybe more careful than I am.
Can one of the Ds maybe defend why they're saying D?
AUDIENCE: [INAUDIBLE]
PROFESSOR: The mean is 100, and when you say-- I mean,
I think that whatever you're thinking is correct,
but I think that the words are a little dangerous.
And why am I concerned about-- you said--
is the mean 100 for all time?
AUDIENCE: [INAUDIBLE] and steady state.
PROFESSOR: And steady state.
Right.
I think that was the-- for long times, the mean number of mRNA
will, indeed, be 100.
So the mean number of m, in this case,
will be Km divided by gamma m, which
is going to be equal to 50 divided by that.
That gets us 100.
Now will it be exactly 100?
No.
It's going to be 100 plus or minus what?
Plus or minus 10.
Right.
Because this distribution at steady state is what?
AUDIENCE: It's Poisson.
PROFESSOR: It's Poisson.
What's the variance of a Poisson distribution?
It's equal to the mean.
So for Poisson, the variance is equal to the mean.
Variance is the square of the standard deviation,
which means that this is going to be plus or minus 10.
That's kind of the typical width of the distribution.
So what it means is that at equilibrium,
we're going to be at 100 and it's going to kind of look
like this.
So this might be 2 sigma, so this could be 20.
But each of these is 10.
So if you want to capture this, you
might want to go out to a few sigma.
So let's say you want to go out to 3 sigma,
then you might want to get out to 130 maybe.
So then, if you want to be more careful
you go out to 140 or 150.
But this thing is going to decay exponentially,
so you don't need to go up TO 1,000,
because the probability's going to be 0 0 0.
Certain once you're at 200 you don't have to worry about it.
But of course, you have to remember the initial condition
we started at 50.
So we started at this point, which means we definitely
have to include that equation.
Otherwise we're in trouble.
Now how much do we have to go to below 50 Any--
AUDIENCE: My guess would be that it would be not much
more than the [? few ?] times 5, because if it were already
at equilibrium that would be the mean.
But it's not, and so the driving force
is still going to push it back to [INAUDIBLE].
PROFESSOR: That's right.
So it's going to be a bias random
walk here, where it's going to be sort of maybe
twice as likely at each step to be moving right as to be moving
left.
That means it could very well go to 49, 48.
But it's not really going to go below 40, say.
Of course you have to quantify these things
if you want to be careful.
But certainly I would say going from, I don't know, 35 to 135
would be fine with me.
You would get full credit on your problem set.
So we'll say-- I'm going to make this up-- from 35 to 135,
134 just so it could be 100 equations.
So I'd say I'd be fine with 100 equations.
So you would simulate the change in the probabilities
of P35 to P134, for example.
So although in principle, the master equation
specifies how the probabilities for an infinite number
of equations are going to change,
you only need to simulate a finite number of them
depending upon the dynamics of your system.
Yes.
Thank you for the question, because it's a very important
practical thing.
AUDIENCE: So in practice, you don't
know what the solution is, which is sort of why you would
[INAUDIBLE].
Do you explain your range and see if the solution changes?
PROFESSOR: So the question is, in this case,
it's a little bit cheating because we already kind of knew
the answer.
We didn't know exactly how the time dependence
was going to go.
How is it that the mean is going to change over time on average?
Exponentially, right?
So on average you will start at 50.
You exponentially relax to 100.
But in many cases, we don't know so much about the system.
And I'd say that what, in general, you can do
is, you have to always specify a finite number of equations.
But then what you can do is, you can put in, like,
reflecting boundary conditions or so on the end,
so you don't allow probability to escape.
But then what you can do is you can run the simulation,
and if you have some reasonable probability to any
of your boundaries, then you know you're in trouble
and you have to extend it from there.
So you can look to say, oh, is it above 10
to the minus 3 or 4, whatever.
And then if it is, then you know you have to go further.
Any other questions about how-- you're actually
going to be doing simulations of this, so these
are relevant questions for you.
All right.
So that's the master equation.
But I'd say the key, key thing to remember
is that it tells you how to calculate
the deterministic evolution of the probability of these states
given some potentially complicated set
of interactions.
Now, a rather orthogonal view to the master equation
is to use the Gillespie algorithm,
or in general, to do direct stochastic simulations
of individual trajectories.
Yeah.
Question before we go.
AUDIENCE: So if we just set it to 0, the probabilities
outside the range we think we need,
would we be losing probability?
PROFESSOR: So the question is whether we're somehow
losing probability.
So what I was proposing before is that you always
want probabilities to sum to 1.
Otherwise it's not our probability
and the mathematicians get upset.
And the key thing there is that you
want to start with-- you have to include all the states that
have probability at the beginning.
So in that sense, you're given an initial distribution,
and you have to include all those states.
Otherwise you're definitely going to do something funny.
You start out with a normalized probability distribution.
And then I guess what I was proposing
is that you have a finite number of equations,
but you don't let the probability leave
or come in from those ends.
And if you do that, then you will always
have a normalized probability distribution.
Of course, at the ends you've kind of
violated the actual equations, and that's
why you have to make sure that you
don't have significant probability at any
of your boundaries.
Does that answer?
Not quite?
AUDIENCE: Because I'm wondering if [INAUDIBLE].
PROFESSOR: So I was not suggesting that you set
the probabilities equal to 0.
I was suggesting that you do what's kind of like what
the equations actually here, which is that you don't
allow any probability to leave.
There's no probably flux on this edge.
So for example, out at P134, I would just say, OK, well,
here's the probability that you have 134 mRNA.
And in principle there are these two arrows,
but you can just get rid of them.
So now any probability that enters here can only come back.
And I've somehow violated my equations.
But if P134 is essentially 0, then it doesn't matter.
So instead of looking at these probabilities evolve kind of as
a whole, we can instead look at individual trajectories, right?
So the idea here is that if we start with the situation--
actually, we can take this thing here.
So we know that at steady state it's going to be 100.
Starts out at 50.
And in this case, with the master equation you say,
OK, well, you start out with all the probability right here.
So you have kind of a delta function at 50.
But then what happens is this thing kind of evolves,
and over time this thing kind of spreads
until you have something that looks
like this, where you have a Poisson distribution centered
around 100.
And this Poisson distribution's going
to be very close to a Gaussian, because you
have a significant number.
So the master equation tells you how this probability
distribution evolves.
Now this is the number m and this is kind of the frequency
that you observe it.
So we can also kind of flip things
so we instead plot the number m on the y-axis.
And we already said the deterministic equations
will look like this.
And the characteristic time scale for this is what?
1 over mm, right?
So this thing relaxes to the equilibrium,
time scale determined by the degradation time of the mRNA.
So these are things that should be really--
you want to be kind of drilled into your head,
and I'm trying to drill, so you'll hear them again
and again.
Now the master equation, indeed, since everything's linear here,
the expectation value over the probability distributions
actually does behave like this.
So the mean of the distributions as a function of time
look like that.
And in some ways, if we were to plot this,
we would say, OK, well, first of all it's all here.
Then it kind of looks like this.
So this is somehow how those probability distributions are
kind of expanding over time.
Now for an individual trajectories,
if we run a bunch of stochastic simulations,
we'll get something that on average looks like this,
but it might look like this.
A different one might look like this, and so on,
although they shouldn't converge there
because that's not consistent.
And if you did a histogram at all those different times
of the individual stochastic trajectories,
you should recover the probability distribution
that you got for the master equation.
So this is a powerful way just to make sure
that, for example, your simulations are working,
that you can check to make sure that everything behaves
in a consistent way.
Now there's a major question, though,
of how is it that you should generate
these stochastic trajectories?
And the sort of most straightforward thing to do
is to just divide up time into a bunch of little delta t's,
and just ask whether anything happened.
So let me--
So what we want to do is we want to imagine we have maybe
m chemical species.
So now these are different m's and n's.
Be careful.
m chemical species, they could be anything, could be proteins,
they could be small molecules, something.
And there are n possible reactions.
And indeed, in some cases people want
to study the stochastic dynamics of large networks.
So you could have 50 chemical species
and 300 different reactions.
So this could be rather complicated.
And these m chemical species have, we'll say, numbers
or if you'd like, in some cases it could be concentrations, Xi,
so then the whole thing can be described as some vector X.
And the question is, how should we assimilate this?
The so-called, what we often call the naive protocol--
and this is indeed what I did in graduate school
because nobody told me that I wasn't supposed to do it--
is that you divide time into little time segments delta t.
Small delta t.
And you just do this over and over.
And for each delta t you ask, did anything happen?
If it did, then you update.
If not, you keep on going.
Now the problem with this approach-- well,
what is the problem with this approach?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Yeah.
Time is continuous.
So one problem is that, well, you don't like discrete time.
That's understandable.
But I'm going to say, well, you know, the details-- a delta t
may be small, so you won't notice.
I'm saying, if I said delta t being small,
then I'm going to claim that you're not going to notice that
I've--
AUDIENCE: [INAUDIBLE]
PROFESSOR: But then the simulation is slow, right?
So there's a fundamental trade-off here.
And in particular, the problem with this protocol
is that for it to behave reasonably,
delta t has to be very small.
And what do I mean by very small, though?
AUDIENCE: [INAUDIBLE]
PROFESSOR: That's right.
For this to work, delta t has to be
such that unlikely for anything to happen.
But this is already a problem, because that
means that we're doing a lot of simulations,
and then just nothing's happening.
How do we figure out what that probability is?
So in particular, we can ask about-- well,
given possible reactions, we'll say with rates rs of i.
So the probability that the i'th reaction occurs is equal to r i
times delta t for small delta t, because each of these reactions
will occur kind of at a rate-- they're going to be exponential
distributions of the times for them to occur.
This is a Poisson process because it's random.
Now what we want to know is the probability that nothing
is going to happen because that's how we're
going to have set delta t.
Well, what we can imagine is, then
we say, well, what's the probability that is, say,
not reaction 1 and not 2 and dot dot dot.
OK.
Well, and this is in some time delta t.
Well, actually, we know that if the fundamental process just
looks like this, then we're going
to get exponential distributions for each of those.
So we end up with e to the r1, and indeed,
once we write an exponential, we don't have to write delta t.
This is just some time t.
For this to be true requires a delta t is very small.
But if we want to just ask, what's
the probability that reaction 1 has not happened in some time
t, this actually is, indeed, precisely equal
to e to the r1t.
Yeah, details.
And this is e to the minus r2t dot dot dot minus.
And we go up to n, r to the nt, because each of those chemical
reactions are going to be exponentially distributed
in terms of how long you have to wait for them to happen.
And what's neat about this is that this means
that if you just ask about the probability
distribution for all of them combined by saying that none
of them have happened, this is actually
just equal to the exponent of minus--
now we might pull the t out and we just sum over ri.
So this is actually, somehow, a little bit surprising,
which is that each of those chemical reactions
occur, and they're occurring at different rates.
Some of them might be fast, some of them might be slow.
The ri's can be different by orders of magnitude.
But still, over these hundreds of chemical reactions,
if the only thing you want to know
is, oh, what's the probability that none of them
have happened, that is also going
to end up-- that's going to decay exponentially.
And this actually tells us something very interesting,
which is that if we want to know the distribution of times
for the first thing to happen, that's
also going to be exponentially distributed.
And it's just exponentially distributed
with a rate that is given by the sum of these rates.
Now that's the basic insight behind
this GIllespie algorithm, where instead of dividing things up
into a bunch of little times delta t, instead what you do
is you ask, how long am I going to have to wait
before the first thing happens?
And you just sample from an exponential with this rate
r that is the sum of the rates.
Maybe it's even worth saying that, OK,
so there's the naive algorithm where you just
divide a bunch of delta t's, you just take a little steps,
you say, OK, nothing, nothing, nothing, nothing, and then
eventually something happens, and then
you update, you keep on going.
There's the somewhat less naive algorithm, which is exact,
so it's not the same concerns, the j hat which
is that you could just sample from n different exponentials,
each with their own rates, and then
just take the minimum of them and say, OK,
that's the that happened first, and then update from that.
And that's an exact algorithm.
But the problem is that you have to sample from possibly
many different exponentials.
And that's not a disaster, but again, it's
computationally slow.
So the Gillespie algorithm removes the requirement
to from those n exponentials, because instead what you
do is you just say, the numbers, or the concentrations,
give all of the ri, give you all the rates.
And then what you do is you sample
from an exponential with rate r, which
is the sum over all the ri.
That tells you, when is the first reaction going to occur.
And then what you do is you ask, well, which reaction did occur?
Because you actually don't know that yet.
And there, it's just the probabilities of each of them.
So the probabilities Pi is just going
to be the ri divided by the sum over the ri, so this big R.
So it may be that you had 300 possible chemical reactions,
but you only have to do two things here.
And they're both kind of simple, right?
You sample from one exponential, gives you
how long you had to wait for something to happen.
And then you just sample from another simple probability
thing here that just tells you which of the n
possible chemical reactions was it that actually occurred.
And of course, the chemical reactions
that were occurring at a faster rate
have a higher probability of being chosen.
So this actually is an exact procedure in the sense
that there's no digitization of time or anything of the sort.
So this actually is computationally efficient
and is exact, assuming that your description of the chemical
reactions was accurate to begin with.
So then what we do is we update time.
This is in some ways-- when you do computations,
when you actually do simulations-- this is maybe
the annoying part about the Gillespie algorithm, which
is that now your times are not equally spaced,
and so then you just have to make
sure you remember that, you don't plot something that's
incorrect.
Because your times are going to hop at different time
intervals.
But that's doable.
You have to update your time and you
have to update your abundances.
And then what you do is repeat.
I think the notes kind of allude to this Gillespie algorithm
but are not quite explicit about what you actually
do to go through this process.
For the simulations that you're going to do in this class,
I would say that you don't get the full benefits
of the Gillespie in the sense that you're not
going to be simulating hundreds of differential equations
with hundreds of different things.
But it's in those complicated models
that you really have to do this kind of Gillespie approach,
as compared to even this somewhat better model, which is
you sample from the different exponentials.
Are there any questions about why this might work,
why you might want to do it?
Yes.
AUDIENCE: What do you mean by sample the exponentials?
PROFESSOR: Right.
What I mean is that you go to Matlab and you say,
random-- I'm sort of serious, but-- sorry,
I'm trying to get a new-- All right.
So you the exponential.
So it's a probability distribution.
So this is the probability is a function of time and then t.
And it's going to look something like this.
This thing is going to be some-- given that, in general, it's
going to be the probability t is going to be e to the minus rt.
And then do I put r here or do I put 1 over r?
AUDIENCE: [INAUDIBLE]
PROFESSOR: Is it 1 over r?
Well, what should be the units of a probability distribution?
1 over time, in this case.
It's 1 over whatever's on this x-axis,
because if you want to get the actual, honest to goodness
probability-- so if you want the probability that t is,
say, between t1 and t1 plus delta t.
If you want an actual probability,
then this thing is equal to the probability density at t1,
in this case, times delta t.
So that means thing has to have a 1 over time,
and that gives us r here.
So this is probability density, and what
I'm saying is that when I say sample from this probability
distribution, what it means is that it's like rolling a die,
but that it's a biased die because it's
continuous thing over the time.
But just like when you have a six-sided die and I say, OK,
sample from the die, you're playing Monopoly,
you throw the die and you get 1, 2, 3, 4, 5, 6.
And you do that over and over again.
Same thing here.
You kind of roll the die and see what happens.
And indeed, you're going to get some practice with probability
distributions on the homework that you're
doing right now because you're asked to demonstrate that you
can sample from a uniform distribution, which something
that's just equally probable across the unit line,
and do a transformation and get an exponential distribution.
And it used to be that everybody knew all these tricks
because you had to kind of know them
in order to do computation.
But now, Matlab, or whatever program you use,
they know all the tricks, so you just
ask it to sample from an exponential with this property
and it does it for you.
But you still need to know what it's doing.
So just to be clear, what is the most likely time
that you're going to get out from the exponential?
0.
It has a peak here but the mean is over here.
Any other questions about how the Gillespie algorithm works?
Can somebody tell me how a protein burst arises?
So we had this original question about whether there
were protein bursts in that model
that I wrote down, where we just had m dot is equal to--
Now what we said was that the master equation would not--
the protein burst would somehow be
there are but you would never see them,
or somehow the protein burst would influence how the mean
and everything have evolved, but you wouldn't actually
see any big jumps.
But then we said, oh, but if you did a stochastic simulation,
you would.
So the claim here is that the Gillespie algorithm,
what I've just told you here, will lead to protein bursts.
When I make that statement, what is it that I actually mean?
If we do a Gillespie of this, will the-- OK,
let's just hold on.
Let me do a quick vote.
Will we have cases where delta n is greater than 1?
If I go through this process, if I'm using the Gillespie
and I'm tracking how mRNA and protein number are changing
over time, will I get these things, protein bursts,
where delta n is larger than 1 in one of these time cycles?
Ready?
3, 2, 1.
So most of the group is saying that it's going to be no.
But again, it's mixed.
So can somebody say why we don't get--
AUDIENCE: [INAUDIBLE] It seems like the structure
of the simulation is to make sure [INAUDIBLE].
PROFESSOR: That's right.
Yeah.
So the simulation as written-- you
could imagine some sort of phenomenological version
of this where you allowed, actually, for protein bursts.
But as kind of specified is that we ask,
what's the time for one thing to happen?
But the claim somehow is, OK, well, we can still
get protein bursts from this.
And how does that happen?
AUDIENCE: You can have the rate for something happening
increase suddenly, and that would happen if we go from m
equals 0 to m equals 1--
PROFESSOR: Yeah, for example, if we didn't have an mRNA before
and we got an mRNA.
What it means that if you look at n
as a function of time during one of these protein
bursts-- before, I was drawing it just hopping up, but really,
in the context of the Gillespie, it
would be that it would hop, hop.
So there would be little time jumps.
So this is a protein burst, but it's really
before this mRNA is degraded, you get 1, 1, 1, 1.
So each of these as is delta n of 1.
So this is whatever, 6, 7.
And then what can happen is that we get the mRNA degraded.
And so then we're going to get a slower thing
where it-- looks like that.
So the Gillespie, everything is being created and destroyed
in units of 1.
But it could be that the time interval over this burst
is just very short, so then it goes up very quickly,
but then it's slower to go away.
So what I want to do in just the last 15 minutes
is talk a bit about the Fokker-Planck approximation.
I would say that all these different approaches are
useful to varying degrees in terms of actually doing
simulations, doing analytic calculations,
getting intuition.
And the Fokker-Planck approach, I'd say it's
more or less useful for different people
depending on what you're doing.
So the basic idea, as kind of you
answered in the pre-class reading,
is that in cases where n is large enough
that you don't feel like you need
to take into account the discrete nature
of the molecules, yet at the same time
it's not so large that you can totally
ignore the fluctuations, then the Fokker-Planck approach
is nice because it allows you to get some sense of what's
going on without all of the crazy details of, for example,
the master equation.
And then it also, because of this idea
of an effective potential, it allows
you to bring all the intuition from that into your study
of these gene circuits.
Now I'm not going to go through the whole derivation,
but if you have questions about that,
please come up after class and I'm
happy to go through it with you, because it's sort of fun.
But the notes do go over it.
I think that's what's perhaps useful to just remind ourselves
of is how it maybe leads to a Gaussian with some width
depending upon the shapes of the production degradation curves.
So the basic notion here is that, depending
on the f's and g's, the production degradation terms,
we get different shaped effective potentials.
So in general we have something that
looks like-- we have some n dot, there's some fn,
and then there's a minus gn.
So for example, for something that is just simple expression,
in the case of-- let's just imagine now
that there is-- if you want we can say it's a protein where
it's just some k minus gamma n.
Or if you'd like, we could say, oh, this is mRNA number.
But something that's just simple production, and then
first order degradation.
The question is, how do we go about understanding
this in the context of the Fokker-Planck approximation?
And it turns out that you can write it
in what is essentially a diffusion equation where
you have some probability flux that's moving around.
And within that realm, you can write
that the probability distribution of the number
is going to be something that-- so there's
going to be some constant.
There's f plus g.
And these are both functions of n.
And then you have e to the minus [INAUDIBLE]
So the idea here is that this behaves
as some effective potential.
Of course, it's not quite true because f and g also
are functions of n, they're are not in here.
But this is the dominant term because it's
in the exponential.
And here phi n is defined as the following.
So it's minus this integral over n of the f minus g
and f plus g dn that we integrate over n prime.
And we're going to kind of go through what
some of these different f's and g's might
look like to try to get a sense of why this happened.
It is worth mentioning that you can
do this for any f and g when it's just in one dimension,
so you just have n.
Once you have it in two dimensions,
so once you actually have mRNA and protein,
for example, you're not guaranteed
to to be able to write it as an effective potential.
Although I guess if you're willing to invoke a vector
potential, then maybe you can.
But in terms of just a simple potential,
then you can do it one dimension, but not necessarily
in more.
And I think that, in general, our intuition is not as useful
when you have the equivalent of magnetic fields and so forth
here anyway.
What I want to do is just try to understand
why this thing looks the way it does for this simple regulation
case.
And then we're going to ask if we change one thing or another,
how does it affect the resulting variance.
So for unregulated expression, such as here,
if we look at the production and degradation as a function of n,
fn is just some constant k, whereas gn is
a line that goes up as gamma n.
Now in this situation, if you do this integral-- and really,
what you can imagine is what this integral looks like right
around that steady state, because that's
kind of what we want to know, if we want to something about,
for example, the width of a distribution.
Well, there's going to b e two terms.
In the numerator there's an f minus g.
In the denominator there's an f plus g.
Now f minus g is actually equal to 0 right
at that steady state, and that's why
it's a steady state, because production and degradation are
equal.
Now as you go away from that location, what you're
doing is you're integrating the difference between the f
and the g.
And you can see that around here these things are
separating kind of-- well, everything's a line here.
And indeed, even if f and g were not linear,
close to that steady state they would be linear.
What we can see is that as you're integrating,
you're integrating across something
that is growing linearly.
That's what gives you a quadratic.
And that's why this effect of potential
ends up behaving as if you're in a quadratic trap.
Now I encourage you to go ahead and do that integral
at some point.
I was planning on doing it for you today,
but we are running out of time.
Once again, I'm happy to do it, just after class.
And indeed, what you can see is that because you're integrating
across here, you end up getting a quadratic increase
in the effective potential.
And if you look at what the variance of that thing is,
you indeed find that the variance is equal to the mean
here.
So what I want to ask in terms of trying
to get intuition is, what happens
if we pull these curves down?
So in particular, let's imagine that we
have a situation where-- I'm going
to re-parameterize things, so again, we're
kind of keeping the number of the equilibrium constant.
But now what I'm going to do is I'm
going to have an fn that looks like this,
and gn looks like-- now gn is going to be some 1/2 of lambda,
and this fn is equal to k minus 1/2 of gamma n.
Now the question is, in this situation,
what will be the variance over the mean?
Well, first of all, the variance over the mean
here was equal to what?
Although should we do vote?
Here are going to be some options.
Question is variance over the mean in this situation.
I'm worried that this is not going to work,
but let's just see where are.
Ready, 3, 2, 1.
All right.
So I'd say that at least broadly, people
are agreeing that the variance over the mean
here is equal to 1.
And again, this is the situation that we've
analyzed many times, which is that in this situation
we get a poisson, where the poisson only
has one free parameter, and that parameter specifies
both the mean and the variance.
So for a poisson, the variance of the mean
is indeed equal to 1.
So the Fokker-Planck approximation actually
accurately recapitulates that.
Now the question is, what will the variance over the mean
be in the situation that I've just drawn here?
So I'm going to give you a minute
to try to think about what this means.
And there are multiple ways of figuring it out.
You can look at, maybe, the integral.
You can think about the biological intuition
to make at least a guess of of what it should do.
The question is, if the production
rate and the degradation rate look like this,
what does that mean for the variance over the mean?
So I'll give you a minute to kind of play with it.
Why don't we go ahead and vote, just
so I can get a sense of where we are?
And also, it's OK if you can't actually figure this out
or you're confused.
But go ahead and make your best guess anyways,
because it's also useful if you can guess
kind of the direction it'll go, even if you can't figure out
its magnitude.
So let's vote.
Ready, 3, 2, 1.
OK.
So it's a mixture now, I'd say, of A, B, C, Ds.
Yeah, I think this is, I think, hard and confusing.
I maybe won't have-- all right.
I'll maybe say something.
It may be that talking to each other won't help that much.
OK, so in this case, what's relevant
is both the f minus g and the f plus g.
And it turns out that f minus g actually behaves the same way,
because at the fixed point, or at the equilibrium,
it starts at 0 and then it actually grows in the same way
as you go away from it.
The difference is the f plus g, where that's very much not
equal to 0.
And f plus g at the equilibrium, this f plus g
here is around 2k, whereas f plus g over here is around 1k.
What that means is that in both cases
you have a quadratic potential.
But here the quadratic potential actually ends up being steeper.
So if this were unregulated, then over here
we still get a quadratic, but it's with steeper walls.
So actually here, this, the variance over the mean,
ends up being 1/2.
It's useful to go ahead and just play with these equations
to see why that happens.
And I think that's a nice way to think about this
is, in this limit, where we pull this crossing point all
the way down to 0, now we have something
that looks kind of like this.
So very, very low rate of degradation.
But then also the production rate
essentially goes to 0 when we're at this point.
So we could still parameterize as k over gamma
if we want, with some-- but we could just
think about this as being at 100 of these mRNAs, say.
But then we're changing the production degradation rate.
And the variance over the mean here--
does anybody have a guess of where that goes?
In this case it actually goes to 0.
And this is an interesting situation, because really,
in the limit where there's no degradation,
and it's all at the production side, what
it's saying is that you produce, you produce, you produce,
until you get to this number, which might be 100,
and then you simply stop doing anything.
You're not degrading, you're not producing.
In that case all the cells will have exactly 100, maybe, mRNA.
And what the Fokker-Planck kind of formalism
tells you is that just because production and degradation
rates are equal, f minus g is equal to 0,
doesn't mean that-- that tells you
that that's the equilibrium, but it doesn't tell you
how much spread there's going to be around the equilibrium.
If f and g are each larger, that leads to a larger spread
because there's more randomness, whereas here, f and g are
both essentially 0 at that point.
What that means is that you kind of just
pile up right at that precise value.
We are out of time, so I think we should quit.
But I am available for the next half hour
if anybody has any questions.
Thanks.