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. 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. 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.
