Flash and JavaScript are required for this feature.
Download the video from iTunes U or the Internet Archive.
Description: Students continued to learn more about Monte Carlo method and further learned about stochastic chemical kinetics in this lecture.
Instructor: William Green
Session 33: Monte Carlo Met...
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 to view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu
WILLIAM GREEN: All right, so today we're going to keep on working on stochastic. And I guess maybe we start-- if you guys have questions from the homework set, since I'm sure you're thinking about that a lot right now. Are there any questions? Should I take this to mean that you have it completely under control? That's good.
All right-- so since we have Monte Carlo integration metropolis totally under control, let's talk about time dependent probability distributions. So in many situations, we'd like to predict what's going to happen in an experiment. And the fundamental equations of that experiment are time dependent, so we have different things happening at different times.
And so we really want to know what d dt of the probability of some observable-- this is going to equal something. Let me fix the lights. All right, so we want some equation like this. And then when you actually made your measurement, you'd be sampling from P.
So probability of the observables-- and so if you made the measurement, you'd have a certain probability that you observe something. And you made a different measurement, you'd have same probability if you made the measurement at this time. So this is going to be probability of the observable depending on the time. And you made the measurements at different times.
So let's conceptually think about this. Suppose we had a box. And it had a uranium atom in it that's radioactive. And it can decay.
So we look in the box. And at time zero, I look in there. And there's the uranium atom sitting there nicely. And I go away for a while.
They come back. And I open the box. And I look again later. And there's some probability that the uranium atom is still there. And there's some probability that it fissioned during the time and it's gone, right?
OK so there must be some kind of equation like this, how the probability that the uranium atom is still there in the box. Right, is this OK-- yeah? So just for that case, what do we think the probability-- how does that probably behave-- any ideas?
So uranium has a half-life or something, right? You guys heard this before, right? So let's just try to think of how this should scale. Whatever it is, it probably should have something to do with a time constant for how fast uranium decays, yeah?
You guys are looking at me, like, so blankly. This is not that hard. It's just one uranium atom. It's OK.
All right, so let's just write a plausible equation. Maybe it's dP dt, where this is the probability that the uranium is in the box. This might be the probability divided by two. Something like that-- would that be-- all right, I guess we need a negative sign maybe.
So then the probability-- the solution to this would be that the probability that uranium is still there is going to be equal to some initial probability e to the negative time over tau. Does that look right? Does this look very reasonable? You might have seen this before, yeah?
Maybe you didn't think of it probability. Maybe you thought about it as number of uraniums in the box. If I had written it as the number of uraniums is equal to the original number of uraniums e to the negative t over tau, you guys would have believed that right away, right?
OK, so I guess maybe this is highlighting a issue here, is that we use a macroscopic thinking. We're used to thinking of things as continuum, right? But if I only have one uranium there, it's either there or it's not there. It's 1 or 0.
So the problem with this equation-- if you compute this, it computes irrational numbers, fractional numbers as the number of uraniums. But of course, the number of uraniums, it's just either-- it's integers, right? There's one uranium atom, or there's two, or there's three. So this equation is not right, right? So this is incorrect.
So that was what you learned in high school, right-- that equation? But that's not right. So this is really the correct equation, that the probability, that you have some number of uraniums that does something. You know-- well this works for one atom, anyway.
So that's the equation for one. So this is probability I have one uranium. It's the initial probability that there's one uranium in a box, which is basically-- my case would be one. I know I put one in there.
And then it just decays. The probability decays. But it doesn't mean uranium decays.
So it means when I do-- if I do the experiment once, and here I have the uranium-- the probability of the uranium is a simple exponential decay-- time. I know I had a uranium atom there when I first put it there. And at later times it's gone away, right? 'Cause at longer and longer times it's more and more likely it would have fissioned. They
Now if I actually do this. And make a-- buy a million boxes. And put a million uranium atoms-- one in each one to start with. And then I just keep coming back once in a while and checking, what I'll see, is that initially I had-- here's the number of uraniums. Initially I had a million uranium atoms. And for some time period I still had a million.
And then something happened and one of them went away, because one of the fissioned, turned into lead-- or whatever they turn into. All right, and so now I had a million minus one uraniums. They live for a while. And then maybe another one went away.
And then this time, another went away. And that time it lasted longer-- right? This is what you really expect to observe. Is this OK?
And so this is sampling from the probability distribution that you would use. Does that make sense? So I'm sampling a million times, so I have a million boxes. This is the probability distribution for one box. And so I'm able to figure this out. Is this OK?
All right, so we'd like to write equations like this for more complicated systems. But you'll still going to have the same confusion, that the probability is going to be continuous variable. And it's going to behave in ways that seem perfectly sensible to you.
But then every time you do the experiment, different things will happen in time. Because you have something at the time-varying probability, and each time you sample from the probability distribution-- every time you do the experiment, it's like a particular instance of sampling from that probability distribution.
And you might get that there is 87 uranium atoms left. And you might get that there's 93. And it's not because of your measurement error. It's because the intrinsic system has its own randomness to it, right?
You know, you look at one uranium atom-- some uranium atoms live for a million years. And some uranium atoms might decay in the next second, right? And you have no idea when they're going to decay. All you know is statistically, on the average, uraniums decay with a certain half-life, which is pretty long, right?
All right, is this totally obvious? I can't figure it out. Is this, like, totally obvious, or totally confusing? It's very funny-- with your eyes and your expressions, totally obvious and totally confusing has the same blank expression.
Can someone ask a question. Maybe that would help us to have a-- this is OK? All right, this is the OK.
AUDIENCE: It's OK.
WILLIAM GREEN: This is OK. All right, so it's totally obvious-- good, hard to tell sometimes. All right, so Joe Scott wrote some brilliant notes trying to explain how do this for chemical kinetics, OK? So you should definitely read them. They're in the-- posted under the materials, general section-- Stochastic Chemical Kinetics-- nice 20 page long thing, something like that.
It's definitely worthwhile to read it. He was very nice at writing very clear notes. So I'll try to explain his notes inexpertly, in a short period of time. I strongly recommend you read them.
So the idea is that the chemical kinetic equations you guys have all used are also incorrect, right-- just like this is incorrect. Because this is assuming that something that's really a discrete variable is a continuous variable.
And we do this all the time, right? We know that material is made of atoms. And so when you measure our mass in a certain volume, there only can be certain discrete values. But we always ignore that. And we always treat is a continuous variable, right?
So you have some flow-- some amount of liquid-- and you have some amount of mass. And it could be any number. And we treat dm dt, or things like dm dx, d rho dx-- we write those down as if everything is perfectly continuous variables. And that's because that in a lot of systems we have, we have so many atoms-- and we can't measure that well anyway-- that we couldn't tell if we were missing one or two atoms.
Whether we end up with an extra half an atom in our equations or not makes no difference, because we're not that precise, right? So we don't worry about it. So we use contiuum equations everywhere and-- but in many cases, what we really should be using are continuous equations for probabilities. And then as going to show you how you do it for chemical kinetics in a minute.
I'd say that this is super fundamental, actually. So quantum mechanics says that everything is wave functions, which is actually like the square root of a probability density. And the real equations are really probability densities. And so that's the real fundamental equations. Really, everything is probabilities. And when you make an experiment, you're just sampling from the probability.
And all our continuum equations are like approximations. They're only valid in, like, the large number of atoms limit, OK? I know this takes maybe a change of thought, because you're so used to those continuum equations you start to believe they're really true, right?
You've probably used them a lot. And you start to feel like you love those equations. They're your favorite equations. You know, you have a relationship with them. Some of them you might hate, some you love. But at least, you have a strong connection with them.
And you just find out that actually they're not even related to the real equations. The real equations are actually something different. This can be very annoying, but it's true. And when you get to smaller and smaller systems that have fewer and fewer molecules in them, then you get more and more sensitive to the fact that the regular equations are incorrect.
And this comes up a lot in small systems. So many of you probably will do research that has to do with biology. In biology there's cells. Inside cells there's subcellular structures. Inside those subcellular structures, there's not very many molecules.
And so you often are at the limit where you have, like, one molecule of a certain type inside this mitochondria, or something. And that's it. And so you wouldn't expect the continuum equations to work.
And in fact there's a big field-- many people at research do in here-- called single molecule spectroscopy, where you can measure the motion of one DNA molecule or something like that. And so you really see the individual molecules, and the effects of the individuality of the molecules.
And this is not only true in biology. For a long time-- there's a field called emulsion polymerization. You guys ever hear of this? So like, all the paint on the walls here is made by that process, where they polymerize methyl methacrylate. They want to make the little beads of polymers all have a small dispersion as possible, to make all the molecules the same.
So what they do is they make a colloidal suspension-- I'm going to draw a little picture. So they have the reactor. They make little bubbles of the monomer that are suspended in some solid. And inside each of these bubbles-- because the way they make them, they're all about the same size. So they all have a lot of bubbles, all about the same size.
And then they polymerize. And they polymerize each of them. And they polymerize. And basically, all the material that's there polymerizes into one molecule. Now how do they make that happen-- is that they introduce, say, a free radical. And one free radical gets into one of these bubbles. And it can polymerize the whole thing, because the reactions are a radical plus monomer goes to a polymer-- bigger polymer, right?
So polymer radical size n goes to polymer radical size n plus 1. That's the polymerization reaction. Now normally, if you did this in the bulk, there would be a lot of radicals. And some of these guys would react with each other. And the radical will be gone, because the two radicals would react, make a new chemical bond, and this process would stop. And so this gives a big dispersion in the size of the polymer chains you get.
So to beat this, they do emulsion polymerization. And they make that number of radicals so small that statistically, at any time, any of these droplets either has zero or one of radicals in it. And there's not-- during the time it does the polymerization, there's not enough time for a second radical to come in.
So it's only got one radical, so that one radical is going to eat up every single monomer molecule in the whole thing. And so it's going take a giant polymer, that's a whole size of the emulsion. If you made the drops too big, you get two radicals in there. And you get some of this happening. And then you get a dispersion. So this is emulsion polymerization.
All right, so cells have small little volumes. Colloids and emulsions have small little volumes. They all have-- that the-- some of the molecules-- the important molecules are so-- have such low concentrations-- the volumes are so small, that the concentration times the volume is less than 1, or about 1. So that's for like-- you know, your copies of DNA in you, in your cells.
It's the same kind of thing for radicals-- free radicals in a polymerization system. And similar things could happen with, like, explosions-- different kinds of rare events, that you might have a very small number of some really important molecule. And then whether or not you have one or not is really important, OK.
All right, so we want to model those kinds of systems, we have to worry about the graininess. So we want to write down what the equation is dP dt. So first you've got to think, well what is P? Let's do a case where I have the reaction of A plus B goes to C, OK?
Suppose this is the only reaction that can happen in our system. The only molecules we have in our system are A's, B's, and C's. And so our probability distribution-- our probability that we want to know is how many A's, how many B's, how many C's do we have in our system, OK? So that's the probability we're interested to know.
We could have 5 A's, 2 B's, and 0 C's. We could have 3 A's, 7 B's, 11 C's-- who knows, right-- different numbers like that. If it gets out to be 10 to the 23, then there's no point. You don't have to do this, because you can just go back your continuum equations, right? But if you have 5, or 3, or 2, or 1, or 0, then you really need to worry about the graininess.
So you want to write down an equation. So here's a probability of n A n B n C. And it's going to equal to reactions which consume-- if I'm in this state that has a certain number A's, B's, and C's, and I can react away, so I only have reactions-- for example, I'm in this state. An d plus B could react with each other, and will react away. So I'll have some reaction. I'm going to write a k, but we'll come back to what this really means exactly in a minute.
So this is the k for A plus B going to C. And we think that this is going to scale with the number of A's and the number of B's I've got, right? And this depends on the probability that I have A B C. So if I'm in this state with n A's n B's n C's-- that probability is going to go down, because the A's and B's can react with each other, right?
Now on the other hand, I have plus the same k. And this is going to be n A prime n B prime. If I have some other state that can react so that it makes this state, then I'll get a plus sign, because probability of this state would increase.
So for example, if n A prime is equal to n A plus 1. And n B prime is equal to n B plus 1 and n C prime is equal to n C minus 1, then that state can react by one of these A's and one of these B's combining together to make a C. That will get me to this exact state here. Does this make sense?
All Right, so that's the master equation if I have a system that could only do the irreversible reaction A plus B goes to C. So the only-- I can go away-- because the probability density of this-- the probability that this state can reduce, because that reaction occurs starting from this state. And probability this state can increase, because I started from this other state here-- this special state here. And that would make the state interest.
All right and so I have equation-- this equation over and over again, for every possible value of n A, n B, n C. Is this OK? Now I have to actually-- remember for chemical kinetics, that whenever you have a forward reaction, you also have the reverse reaction. So there's actually a couple more terms.
So now I have another loss term. This is C goes to A plus B. This is going to tell of the number of C's. And that's P of n A n B n C. And then I have another source term. And this is where I start from the state with one more C. Is that OK?
So this set of equations is called the kinetic master equation. And this is the equation that describes the real system. So this is the real probability of finding the system with any number of A's, B's, and C's. And I have an equation like this for every different possible number of A's, number of B's, number of C's. All right?
And it's not bad as an equation goes. It's a linear differential equation, because it's linear in P's, right? There's just first-- everything's first power in P. So in general, I can rewrite this as dP dt is the some matrix times B. And the matrix entries are these prefactors-- things like this and this-- sorry, this, this. These are the different elements in the matrix, m.
So you should go ahead and be able to write any kinetic equation system-- rewrite it in this form if you want. It's called the kinetic master equation. Now if I do an experiment like I talked about with uranium, where I put one uranium atom in there.
And all uranium atoms can react to do is fall apart to lead plus neutrons, or whatever the heck they fall apart to. So I can write down-- there's like one reaction that uranium nuclei do. And so I can write it out. It will just be something like this.
And this equation will work for a case where I put one uranium atom in there to begin with. I could do the case where I put two uranium atoms in there to begin with. I could write this down. Is that all right?
OK, now if I had a case where I had 27 A's, 14 B's, and 33 C's to start with, then you can see that the number of possibilities here might get pretty large. 'Cause I have to consider, well, I could lose one A, and lose one B, and get one more C.
But then I'll have an equation for this term which will look just like this. But it'll go to ones where I have lost two A's and two B's and gained two C's. And then all the way down, until you get down to one of the n's to be zero.
And so it could be quite a lot of possible equations, right? So the dimension-- the length of the vector P might get really long, right? So the length of P is something like the max number of A's, max number of B's, max number of C's.
All right, so if I have-- I'll allow-- consider systems with up to 10 A's, 10 B's, 10 C's-- then this is 1,000-- P's a length of approximately 1,000. And it might be less than that, because maybe there's some of them-- some particular states I can't get to. So maybe it's less than 1,000, but something like 1,000.
On the other hand, if I allow it to a million A's, a million B's, and millions C's, then I, like, 10 to the 18 length of P. So P gets pretty long, pretty fast when the number of possible states goes up. Is that all right?
So-- now, it's linear though. And it's linear. And this is constant coefficients. And so this is really not-- and sparse. So this is really-- it's not that bad. And in fact Professor Barton's group has worked on these kind of problems. And he has to solves that can work with P has dimension up to 200 million.
So if you can keep the number-- say it's 200 million or less, you can actually solve this exactly. Now solving this exactly is really good, because then you actually have the exact probability of any observable measurement that you would possibly make at any time, all computed, OK? So that's really great.
But as I just pointed out, if I make these numbers, like, 10 to the 6, 10 to the 6, and 10 to the 6-- this is like 10 of the 18. That's a lot bigger than 200 million. And so Professor Barton's numerical methods are going to poop out at some point, and not going to work for you anymore.
So there's a limited number of problems you can do. So if your problem small enough, you can just solve this-- for really small ones you can solve this with OD45 or OD15s But because it's linear, there's special solution methods. I don't know if you remember solutions of linear equations this, you can relate them to the eigenvectors and eigenvalues of m. And so you can do-- that's one way to do it.
And then there are special methods you can do. Talk to Professor Barton. If it's sparse-- and it has certain kinds properties, you can make this really-- you can make quite large systems, up to 10 to the 8 kind of size.
If we get a system with more than 10 to 8 possible states, that means more than 10 to the 8 possible things could be observed or the result of experiment could have more than 10 to 8 possibilities, then this is-- you're not going to be able to solve it this way. And so then we're going to have to use a stochastic solid solution method to try to figure-- to do-- so we're going to try a sample from P of t.
So I guess I should comment, again, I wrote it this way. But the thing-- P is definitely-- it depends on time. So it's you'll have different P for every time. Yeah?
AUDIENCE: [INAUDIBLE]
WILLIAM GREEN: Yeah maybe I'll just do better to-- I can substitute it in here. Is that all right?
AUDIENCE: So, like, if you transfer all the neighbors in it, that's--?
WILLIAM GREEN: Yeah, 'cause a key idea is that chemical kinetics only does one step at a time. So if you have a reaction like this, at any instant the probability that two reactions are going to be happening exactly the same instant is negligible. So at any instant you're only changing the numbers of the atoms by one unit, yeah.
AUDIENCE: Then shouldn't you add another A minus 1 and an E minus 1?
WILLIAM GREEN: Maybe I got this backwards. n A plus 1 plus 1 minus 1, minus 1 minus 1 plus 1-- I think got-- I think this is right. So this is for the forward reaction that always reduces the number of A's and B's. And this is the reverse reaction that always increases them.
And so I think you generally get like this. You always get-- for any reaction, you'll have four terms-- so the forward and reverse, starting from initial state, and the forward and reverse that make your state. And then if you add more reactions, you have a system where you've got a lot more reactions. Maybe you have C's that can just decompose into B's or something-- different reactions. Then you'd have a sum like this with terms like this from every reaction. And you'd go through all your list of all your reaction.
I guess there's one funny thing I should point out to you about this, which is very important for this emulsion polymerization case-- Is that if you have a reaction of the type A plus A goes to-- I don't know-- B. That might be a reaction you could have. That when I write down the equation, if I only have one A in my sample, this probability of this happening is 0, right? Because I need two A's for this to happen.
So really, the equation for this case is n A times n A minus 1. So it's k A plus A goes to B times n A times n A minus 1. And this is crucial for, like, the emulsion polymerization case, 'cause the fact that I don't have the second one is really important. Otherwise I always miscompute that this radical could react with itself and terminate itself, which is not true. Is that all right?
I commented that these k's are a little bit odd. So if you remember, if you have a k for a unimolecular reaction, like C goes to A plus B, what units would it have? What are the units for regular kinetics?
AUDIENCE: Inverse seconds.
WILLIAM GREEN: Right, inverse seconds, right? So you would normally say this thing has units of inverse seconds. And that works fine here, right? Because I have dP dt, and it's in for a seconds times P, so it's fine. However, this guy-- to make the units work out, still has to also be inverse seconds, but that's not normally what you would have for a bimolecular reaction.
And the other thing to watch out for here is that I have n's. But normally you would write them with concentrations, right? You have the-- you'd would write something like k times the concentration of C. It's the normal way you write your equations, the right equations.
But this is not the same, right? This is n. This is unitless. Now it turns out, in this case, actually, you can use exactly the same as k and it works. But for the bimolecular case, it definitely does not work. The units are wrong.
And also this concentration thing implies that we think the rate really has-- suppose we have k A k B. I could rewrite this as k times n A over the volume, n B over the volume. Is that all right? Actually one thing I should warn you about-- that usually when I write this, I write moles per liter. So I have factors of 10 to the 23 between this and this, because these n's have got to be the individual molecules for this to work.
So watch out. First of all, there's a 10 to the 23 somewhere in here, right-- Avogadro's number. And then this is what-- how we would normally write it. But this k has the wrong units for us here. And so this guy here is the normal k divided by the volume.
I think that's right. Maybe it's the other way around. No, this is normal k times the volume. That's right. So normally we would write units k normal-- normal kinetics, macro-- k for A plus B reactions, the units are centimeters cubed per mole second.
But we want them to be units of-- the new k has got to be units of per second. And so the new k is equal to the old k new old times the volume. Is that right-- no divided by the volume. I was right the first time. This is volume-- yeah. Thank you-- centimetres cubed, centimeters cubed-- yeah, right. Cancelled out, good for a second. So this is k per volume normal k divided by n-- correct.
And you may wonder, well how come it's only volume to the first power? Before I had volume to the second power in the denominator. And it's because I'm using n's up in the numerator. And it's just like here. I'm using n's instead of n over v. So I've sort of multiplied through by v everywhere. That's one way to look at it.
All right, no questions yet? Yeah?
AUDIENCE: So when you [INAUDIBLE]
WILLIAM GREEN: Yeah, so I need a 10 to the 23 as well. So there's got to be an Avogadro's number. And you'll have to help me where the Avogadro's number goes. But it's got to be in there somewhere, to get from moles to molecules. So the numbers are a lot different. Is this OK?
And this factor, this n minus 1 thing-- again, when these are macroscopic numbers like 10 to the 23, this makes no difference, right? Because it's minus 1. So we can fit our data ignoring the 1-- It's 1. Perfectly well, the law of mass action-- it works fine.
So you never notice this unless you get down to where n is very small. But this is actually the truth. This is really what the correct equation is.
OK so now we have these-- the true equations for predicting how the probability distribution evolves in a reacting system. By the way, you can add transport in as well if you want. So you have control volumes, and you know formulas for the rate at which stuff is transported in and out of the control volume. You can add those terms in as well. And they'll affect the probability distribution.
No problem, you should be able to do the same conversion. So just-- we've converted the reaction equations. You can write a similar thing. Once you get right the macroscopic transport equation, you can figure out how you would reduce it down to a time constant for transport that would work for the-- for this microscopic version of the equations.
So transport reaction-- you can write the whole thing down as a master equation. That's a real solution, the real equation. And the problem is that the real equation oftentimes is too big, because P scales up so rapidly with the number of possible states.
So then we have to figure out how are we going to solve it. And the way to solve it is by-- called kinetic Monte Carlo. And this method was invented by a guy named Joe Gillespie. And the journal paper is posted.
And this is just like with Monte Carlo-- Metropolis Monte Carlo that you just did. We're trying to find a way to sample from the true probability density without ever computing the true probability density. So in the Monte Carlo integration you're doing for homework, you're sampling from the true Boltzmann probability density, but you never actually supposedly write down what that P is, right?
And so we're going to do the same thing here, is we're going to-- from the differential equations-- from this differential equation system here, we're going to try to somehow never ever solve this directly, but instead just sample from the solution, P, which depends on time. And the way it's done is really to try to follow individual trajectories of what could have happened, OK?
So at any time-step, we have some number of A's, B's and C's. And we're going to say, well in the next time period, some small time period, what could happen? And we're going to make the time period so small that the only thing that can happen is nothing happens, and A, B and C stay the same, or one event happens which will change one of the numbers, the A's, B's, and C's, by the same one reaction amount, OK?
And then after we accomplish that, now we have a new state-- A, B, and C with some new values. And we'll do the same thing again. And we just repeat that over and over. So it's very brute force. You're just doing, like, what a single system would do.
So for example, for the uranium atom case, we have one uranium atom in there. We wait for some time period, say 20 minutes. We compute the probability that the uranium would have decayed in the last 20 minutes. It's pretty small.
And then we'll say, OK, now we'll draw a random number. If the random number is less than that probability, then we think that the uranium is still there. And if the random number is bigger than that probability, then the uranium is gone. If it's still there, then we'll do it again. And we'll just keep doing that over and over again.
And we'll sample what could happen. The uranium atom could live for 38 years and 14 minutes, and then disappear. And that would be one example. And then we go back and we do it again from a different trajectory.
And we do it over and over and over again. And we do it over and over again. Those trajectories are sampling from what could happen, with the right probability density. Does this make sense?
So the way Gillespie says to do it is figure out the expected time until something would happen. We'll call this time tau. And we'll pull-- actually let's not use tau, because that's the only one that-- I'll try to be consistent with Joe's notation. He calls this thing 1/a.
So a is like the expected rate per second. And 1/a is the expected time until something would happen. And then we'll draw a random number r-- he calls it r1. So we just call the rand function. It gives me a number between 0 and 1. And I'll say that my time-step, delta t-- which Joe calls tau, is 1/a the logarithm of 1/o.
So I think there should be a negative sign. No, 1/r is positive. That's right-- positive. All right, so I choose a random number. Suppose it's 1/2. Then this fraction is 2-- 1 over 1/2 is 2. Logarithm of 2 is some number bigger than 1.
And I pull that. And I get some-- actually, some number less than 1-- logarithm of 2-- anyway, some positive number. Multiply it by the expected time until something happens-- and that's the time, in this particular simulation, when the next thing happened. And then I repeat this again, and again, and again. And I get the period of time intervals between when something happened.
Now every time something happened, I have to figure out, well, what happened? So if I look at my-- plot back here when I had the uranium, I waited a long-- in this particular case, here's my tau. I'll call this tau 1-- the time that it took for the first thing to happen.
And then, well-- what could happen? The only thing can happen with my single uranium atom case is it decayed. So it went away. All right, and then I do the random thing again. And I draw. And I get a different number this time. How long it is 'till something happened.
Now notice-- maybe it's not so obvious-- that in the time until something happens changes as you run. So for example, if I started out with 100 uranium atoms in 100 little boxes. And I'm watching them. After I've run this for a long time, I only have 50 uranium atoms left.
And that means that the average time period until something happens will be slower, because I don't have as many. The more I have, the more rapidly something will happen. Does that make sense?
So I have to compute this a thing at each iteration. So I compute the time until something happens. Then I'm going to compute what happened. I'll tell you what that the second. In the uranium case, the only thing that can happen is one of the uranium atoms decayed. So I took that away. Now I have to compute a again.
So a is the sum of all the things that can happen. So in my case where I have 100 uranium atoms, I have 100 different boxes. In each one of them, there is a rate at which uranium atom probability is decaying. It's like 1 over the normal half life uranium-- or something like that.
And I add them up. So a is the sum of the rates of everything that can happen. So we just list out all the different things that can happen. So in the case-- suppose I have two boxes-- two uranium atoms. I can have decay of uranium atom number one or uranium atom number two.
So I have-- this thing would be equal to-- so two uranium atoms-- a is equal to 1 over tau normal for a uranium atom plus one over tau normal. Because they're both normal uranium atoms. So it's equal to 2/tau.
And so 1/a is half of the normal lifetime. And that's when I would expect that something would happen by then, OK? One of them would probably decay.
So that would be the a value I would use here. Draw my random number, get my tau until something probably happened, then figure out what happens. Now in the kinetics case, several different things can happen, right?
I wrote down if I'm in state n A n B n C, I-- the state can change by two different ways. I can move out of that state this way or this way. Two different reactions can happen-- a forward or reverse reaction. And they would both take me out of the state. So those are two different things that can happen.
And so I have to add those rates. So I would add this number and this number. And they would be the two things that could happen. And that would give me the a value.
Now I figured out something happened. But now I'm not doing uranium. Now doing, say, number of A's. Now I figured out something happened. But now I have to figure out what happened. And the number of A's could have gone up or gone down, because if it was a forward reaction, the number of A's went down. If it was a reverse reaction, the number of A's went up.
So I have to decide which way is going to happen. So I could have either the number of A's went down, or the number of A's went up. Those are two possibilities. I have to figure out which one happened. I have to figure-- draw another random number to figure out which one probably happened.
So I'll list out the list of all the things that can happen with their rates. And one of them, say, is twice as big as the other. So that'll be twice as likely to happen. So I have to choose a random number, and use that to choose the next one.
So first I figured out-- this is actually what we did first. I figured out what my rates are. Second thing is I figure out the time until something happened. Third thing is I'm going to figure out what happened.
And what I do is I say-- I draw a random number called r2. And if there's only two possibilities like in this case, then I can say if r2 is less than the rate of process one divided by the sum of the rates, which is the same as the a, right? So this is-- it this is true, then process one happened.
On the other hand, if the rate-- if the random number I chose is larger than this ratio, then I'll say process two happened. And that's how I decide whether a number of A's went up or went down. And then I would recompute the a-- depending on what happened.
Suppose I do this. And I found out that the number of A's went up-- because I-- by drawing the random number. Then I'll go and I'll recompute the A. Now the rates will all be different, because n A is bigger. So the rate of the n A reaction is faster. So I recompute it.
I redraw a number here, and get some different time, tau 2. It's not the same as tau 1. And then something happened here. I do the same process here out what happened. Maybe this time A reacted away and it went back down.
And then I do it again. And this time I draw a longer time, tau 3. And I do-- what happened? And maybe this time, again, something reacted away. And it went down.
And on, and on, and on-- and after my computer is done doing that for a while, like, after a long enough time, I say, OK, I'm done. Oh-- I'm not really done. I have to go back and start it again, and run another simulation to run through all the possibilities.
And I want to save all these results of all these simulations, because these guys are samples from the probability distribution. And if I built histograms of stuff from those guys I can figure out what really is going to happen. All right, so it's a way to represent the P. Questions?
All right-- perfectly clear? All right, so what's nice about this is actually, it's not that hard once you do this once. It's not hard to do it. And you can do it for any problem at all. And so it's like-- often cases people do this even for cases where you actually could solve-- you might be able to solve this equation. A lot times, people won't do it. They'll do this instead just because it's so easy.
So though maybe it looks hard to you right now, it's only a three step thing. And just do it over, and over again. You're not doing it. The computer's doing it-- no big deal.
And so this one you have to, like, look at the sparsity pattern of your m-- and some complicated stuff. So this is actually done extremely commonly. All right, see you on Friday.