Stochastic Modeling

• Download this transcript - PDF (English - US)

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.