PHYS-413 Topics in the Physical Basis of Physiology

From McGill University Physics Department Technical Services Wiki
(Redirected from PHYS-413)
Jump to: navigation, search

Hi everyone.

Some people are having trouble logging in. If so, email me the username and password you would like to use, and you will be added to the database. - Adriano ()

For now, you can edit this page on the actual Wikipedia (takes 15 sec to make an account, though it's not even necessary), and I will transfer the contents here with your name attached.

Also, note that anything added is considered a contribution, even a one-line question. If you think you know something (even just a vague general idea), just post it, and we'll see how much we can get out of it.



Homework 3

I am having a little trouble with 1.(vi) (isn't the noise 0, since we are approximating dm/dt =0?) and 1.(viii). Also, even after correcting the code typos, none of the variance graphs make sense. For (vi), you can use the result of equating coefficients for the x terms, and you see that it is almost exactly the same form as the equation used in the notes. The ODE can be solved for F(x), and you can find and in terms of and .

Q.1 Stochastic Gene Expression

  • A few minor sticking points I had. As usual, easy in retrospect, but so it goes.
    • and but
    • for vii) we can arrive at the expression by setting , and then use to approximate the term as . Still figuring out the reasoning behind this though...
  • I am still stuck with (viii)... The equation for that I had from (iv) is:
    • .
  • I tried solving it for , and plugging it into , but I can't get that expression. Any ideas, or is this to remain a monologue?

What equation are you talking about? The one that we arrived at from equating coefficents is this: , and can't see a way to solve this into the form of the answer.

Q.2 Langevin Equations

  • (i) Steasy-state:
  • (i) Langevin Eqns.:
    • (or simply the mass-action equations + the xi's)

"Excercise 2 (i): I wonder what happened to the v0 term in the x2 - equation...?"

I am not sure, I'll check later. By the way, I think the way I did the Langevin eqns was excessively complicated. It might be enough to just write down the equations from mass-action, and add the corresponding random variables.

this is my first equation

and the random variables should be accorded in order:

  • , ,

okay I kinda get that you have to do the taylor expansion on steady-state, but how? can you explain why f1*A*M becomes f1*(Msx1+Asx2) ?


  • (ii) a) I found the following non-zero rates (z,m,a):
  • The cross-correlations I get are:
  • One thing I am pretty sure about is . Note, I am also getting . not sure which is right...
  • I am pretty sure the signs are right at least. I'd like to check these against someone though... anyone? I can add details for this, if needed.

Q.3 Simulating stochastic gene expression

I am getting very weird results for the variance: variance of M vs. Tr, for 10 runs. VarMvsTRruns10.gif <--- What the hell is this??

Ok, so I am not the only one: "There's obviously something wrong. I tried it myself and I got something weird too. The mean seems ok, but I can't figure out what is going on with the variance. It looks like xpp has some trouble computing it... I'll think about it. For the moment, I guess you can just "guess" by looking if variance increases or decreases depending what you're looking at. Catherine"

Hey, some of us produced an ode file for 3. v) doesn't seem to work. Any ideas? Sorry about the % for comments...# just numbers the lines, little confusing (fixed it. putting a space before the line leads to the 'boxed code' look). There is definitely a problem with our f1 value. We are trying to get the units to agree.

Did you guys try setting f1=b1=0 to check if the rest of the code is ok (never mind, I just did and it looks fine)? I'm also trying to figure out the f1 units...

Where did you guys get the value f1=0.0083 from? The 0.0083 came from Emma, but we are not sure where it came from. We are trying to put into units of 1/s, so messing around with Avogadro's number and the cell volume. No success so far...I suspect the code is quite busted, as Catherine said. Alright, so the 0.0083 comes from the following: f1 * particles/(Avogradro's Number * Volume) = particles * (1*10^7)/((6.022*10^23) * (2*10^(-15))) = 0.0083 * particles

When I talked to Catherine at some point in the past she seemed to think the code was fine and that the variance problems were coming from the way that xpp was handling memory issues...So I suspect that we're just getting screwed over here by an unclear explanation of where the parameter values come from, not pre-broken code. I'll tell you if I come with anything...

Alright, so I think I've found the problem with the code you have (I initially made exactly the same changes to the code as you did and also got craziness, but now at least I'm getting something that integrates properly). For some reason Catherine wrote into her code for defining variables like M' that they should be the maximum between either 1 or a given formula. The thing that is weird about this is that the variables M', N', etc. define the values for M, N, etc. for the next step, so her code just artificially makes sure that the value of mRNA or protein never falls to 0 particles in the following step. I'm not sure why she did this, especially since you can take out the "max" part and the simulation works fine - either way, what matters is that you can't use the max statement with the values of D and C once you have the negative transcriptional feedback system in place - because you need to always have C + D = 1 (i.e. DNA is either free to be transcribed or is bound up in the complex), and using the "max" statement allows C + D > 1 (you can check the values of C and D in your code by runnning a simulation and you'll see that they immediately sum to 2 in the first time step, and can then just get higher). So the solution is to just use the following lines of code instead:

% D'=D-z5+z6
% C'=C+z5-z6

Unfortunately this gives a result that is not as nice-looking as what we had before. But at least it's physically valid...

Kick ass. It works nicely. Yeah, DNA is always in D or C, so the sum has to be conserved. Sameer has something he would like to say to you:

Jeff (corrected from Adriano), I want to have your babies. -Sameer.

Alright, this is highly immodest of myself, but I was the one who posted all the stuff about the code corrections - Jeff

P.S. I'm clarifying this point because I don't want Adriano to get all the Sameer baby-making action to himself

% gillesp_gene.ode
% gillespie algorithm for gene expression
% D -> M + D (v0)
% M -> M + N (v1)
% M -> (d0)
% N -> (d1)
% N + D -> C (f1)
% C -> N + D (b1)
par d0=0.00385, v0=0.01, d1=0.0001925, v1=0.04, f1=0.0083, b1=0.1
init M=0, N=0, D=1, C=0
%  compute the cumulative reactions
%The units of DN were changed to Molarity from number of molecules because of the f1 units
% choose random number between 0 and p6(=a0)
% Decide which reaction occurs next
% zi = 1 if reaction i occurs. If s2<p1, reaction 1 occurs
% If p1<= s2 <p2, reaction 2 occurs and so on.
% Time for the next reaction
% Update the molecules number depending on which reaction occured
@ bound=100000000,meth=discrete,total=72000,njmp=1
@ xp=tr,yp=M
@ xlo=0,ylo=0,xhi=72000,yhi=10

1. Biochemical Reactions

1.1 Law of Mass Action

The "law" is pretty straightforward. What I am confused about is the difference between the "equilibrium approximation" and the "quasi-steady state approximation". Any one have an explaination?

1.2 Enzyme Kinetics

The simplest reaction scheme is:

1.2.1 The Equilibrium Approximation

Instantaneous equilibrium between the subtrate and complex implies (??) . Also, we have , so finally we have:

where . The resulting velocity of the reaction is:


  • Linear at small s
  • Saturates to at large s

1.2.2 The Quasi-Steady-State Approximation

  1. First, introduce dimensionless variables (keep an eye out for any small parameters, such as in this case).
  2. The quasi-steady-state approximation is that .
"Not the same as taking , but is equivalent to ." Which confuses me at this point :/

1.3 Glycolysis and Glycolitic Oscillations

The actual process seems to be a complete mess. But as far as I understand it, we mainly focused on the mathematical analysis of the resulting equations.

Sel'kov Model

After non-dimensionalization, we find that and are both "fast" variables, since (?) they are of the form This means they can be set to their quasi-steady-state values. The result is:

Now the idea is to show that this has oscillatory solutions. The first step would be to find the steady state solution. Then we'd linearize the ODEs about the steady-state, and the we examine the eigenvalues of the resulting linearized system.

  1. Find the nullclines (), then find their intersection.
  2. Linearize the equations (Taylor expand f about fixed pt., and truncate.). The result is:
where , and denotes the deviation from the steady-state value of . The solution to a linear system is: Plugging this into the equations and simplifying, yields:
which, to have a solution, the determinant of the coefficient matrix must be 0. This gives us the characteristic equation:
Finally, we see that for , we have , and vice versa. At H=0, we have eigenvalues:
which is oscillatory with frequency .

1.4 Math Background

1.4.1 Basic Techniques

Rules of thumb for nondimensionalization:

  1. Start by rescaling the independant (time, space, etc.) and the dependant (concentrations, etc.) variables by "typical" units (usually so they are order 1).
  2. Rewrite the equations using these scaled variables, and find dimensionless combinations of the remaining parameters.
  3. Note which dimensionless numbers are much larger or smaller than the others.

1.4.2 Asymptotic Analysis

1.4.3 Enzyme Kinetics and Singular Perturbation Theory

2. Cellular Homeostasis

2.1 The Cell Membrane

This may be basic to most of you physio people, but it's the first I've heard of this stuff.

Stuff in the cell membrane:

  • Lipid Bilayer
  • Aggregates of globular proteins (free to move within the layer)
  • Water-filled pores (~0.8 nm)
  • Protein-lined pores (channels), which allow passage of specific molecules.

Stuff outside and inside the cell:

  • Both in and out contain dilute aqueous solution of dissolved salts (mostly NaCl, and KCl).

How stuff gets in or out of the cell:

  • Passive transport:
    • Osmosis: the process by which water is transported thru the membrane.
    • Simple Diffusion: diffusion of small molecules thru the pores (e.g. water, urea, hydrated chloride ions, oxygen, CO2, etc...).
    • Carrier-mediated diffusion: a molecule binds to a carrier molecule that moves thru the membrance (e.g. glucose and amino acids are carrier-mediated).
  • Active Transport:
    • Energy required to maintain (usually from ATP→ADP)
    • Some examples: Na+-K+ pump (Na+ out, K+ in); Ca2+ out (using ATP); Na+ in Ca2+ out (using energy from the concentration gradient).

2.2 Diffusion

2.2.1 Fick's Law

where u is the amount of some chemical species, and f is the local production (source term). If D is a constant (usually the case), we have:

2.2.2 Diffusion Coefficients

for sperical solutes, large compared to the solvent. Here k is Boltzmann's, T is temperature, is the coeff. of viscocity for the solute, and a is the radius of the solute molecule.

2.2.3 Diffusion Through a Membrane: Ohm's Law

We have two resevoirs (cl at x=0, cr at x=L), and assume 1D. Therefore:

The steady-state solution is:

From Fick's Law, we have:

which is a constant. The quantity D/L is a conductance per unit area.

2.5 Active Transport

2.6 The Membrane Potential

2.6.1 The Nernst Equilibrium Potential

At equilibrium the potential difference, VS, across the membrane is given by the Nernst potential,

where R is the universal gas constant, T the temp., F is Faraday's const., k is Boltzmann's, q is the unit charge, and z is the charge on the ion S. Note that .

2.6.2 Electrodiffusion: The Goldman-Hodgkin-Katz Equations

The Nernst-Planck equation describes the flow of ions due to both concentration gradients and the electric field:

or, in 1D:

If we solve these equations in the domain , and label the quantities as subscripted i or e (for internal or external quantities) and then integrate these equations, we recover the Nernst potential.

The constant field approximation

If we assume a constant electric field we have , and we can solve the resulting equations. Using the usual boundary conditions for the concentrations, we get the Goldman-Hodgkin-Katz (GHK) equation (for the current):

where is the permeability of the membrane to S.

2.7 Osmosis

2.8 Control of Cell Volume

3. Membrane Ion Channels

3.1 Current-Voltage Relations

3.5 Channel Gating

4. Excitability

4.1 The Hodgkin-Huxley Model

  • Space-clamp technique: eliminates axial variable.
  • Voltage Clamp: step voltage up to some fixed value Vc, and measure the ionic current required to maintain this voltage.

4.1.2 Voltage and Time Dependance of Conductance

  • Separating Ionic currents
    • Seawater ↔ Choline
    • , where j=1 is for seawater, j=2 is for choline.
    • . Also, , and so we have:
    • So, we can solve for and , since and K are determined experimentally.
  • The conductances are therefore (since data supports linear I-V):
which are time-dependant.
  • The Potassium conductance
    • We can assume satisfies , and (through "luck or genius" as Mackey said), assume further that for some variable n which satisfies a first-order ODE.

4.1.3 Quantitative Analysis

The fast phase-plane

  • Fix the slow variables n and h at their resting states, and consider only the behavior of m and v.

The fast-slow phase-plane

  • Assume m is an "instantaneous function of v" (whatever that means) and "thus" at all times.
  • Note that , so we can eliminate h.
  • The Hodgkin-Huxley model is now two-dimensional:

Where the hell is the cubic coming from? Is it an actual cubic dependence, or are they simply noting that it "looks cubic" in order to connect it to the FitzHugh-Nagumo system later? update: I am beginning to suspect strongly that this is the case... how confusing.


Anyway, I finally sort of understand the reason for the trajectories we were shown (thanks to Keener/Sneyd). "Because v is a fast variable and n is a slow one, the solution trajectories are almost horizontal except where ." (i.e. except near the v nullcline). I am not sure if this was mentioned and I just missed it or didn't understand it, but anyway, there it is...

4.2 Two-Variable Models

8. Passive Electrical Flow in Neurons

8.1 The Cable Equation

In the H-H examination of neurons, they used the space-clamp technique to eliminate the space dependence. Now we examine what the space dependence is in neurons.


From Kirchoff's laws, we find that:

  • The axial voltage drops across a single element are:
where and are the (internal and external) resistances per unit length. Also, confusingly, the currents and (internal and external) are actual currents, while (the transmembrane current) is a current per unit length.
  • Adding the currents going into and out of an arbitrary vertex, we have (directions are confusing, but idea is simple):

In the limit , we have:

If we introduce , and use, and rearrange stuff, we get that:

More manipulations yields:

Finally, is the sum of ionic and capacitive currents (and maybe some applied):

so, we have:

Assuming and are constant, and defining a whole bunch of constants yields:

where ri and re are the intra- and extracellular resistance/length, Rm is the membrane resistivity, and p is the axon perimeter.

8.2 Dendridic Conduction

Assuming the membrane is an Ohmic resistor (V=IR), we have:

With some new variables (, ), and assuming the membrane is an Ohmic resistor, we then have the linear cable equation:

We saw in class that, for certain conditions (semi-infinite & voltage clamped ? don't have notes for that day) this has the solution:

where .

Unsorted Topics


  • Problems 1 or 2: Ion transport problems involving membranes (50 pts)
  • Problems 3 or 4: Chemical Kinetics & Cable Equation (50 pts)
  • Problems 5 or 6: Bonus (25 pts)

CLARIFICATION: Is this an open book/notes exam??

What an optimistic question... why did I think he would allow open notes? Anyway, no books, no notes...


membrane stuff :

chemical stuff :

hope this helps for mid-term (-zeshan)

Homework 2

Problem 1

In my lecture in class I derived the cable equation in the form

and then claimed by a change of variables and with , I could rewrite the cable equation in the form

Verify that this procedure is correct.

Problem 1 - Comments

Straightforward... I went backwards though (plugged in the definition of U into the last equation and showed it equals the first)

Problem 2

In class I derived the solution of the cable equation and it involved the use of the ‘error function’.

(a) Verify that the solution that I wrote down on the board is indeed the solution of the cable equation.
For those of you that (like me) didn't have the equation in your notes, here it is (thanks to Zeshan):
Where .
  • I just differentiated twice in X, and once in T and showed that they are equal.
  • Note that . Use this and the chain rule to get the derivs you need.
(b) Sketch or plot as a function of for .
  • I got something like this:

V(1,T).gif V(2,T).gif V(3,T).gif

(c) Transform back to the form , take the limit of the solution, and show that it reduces to the solution of the time independent equation
  • For this part, all you need to know is that

Problem 3

This problem is causing me the most confusion, since I am not 100% sure of the questions. I'll post my thoughts on this soon. Anyone else having problems with this one?

Ok, we seem to making a little progress. Update coming soon. Anyone else have any input?

Wow, (a) and (b) are really easy in retrospect... just multiply both sides of the HH equations by and define the scaled time as (note that this mens that is the derivative with respect to T). From there it's clear which constants are what. (actually, this tau thing is causing me confusion too).

(c) is still causing trouble, because we are not sure what we can assume. As we understand it, from the notes and other sources, one can make the "observation" that v and m are fast variables and "lump them together" into a new (?) variable. Also, h and n are slow variables, so we "lump them" into a variable w. Now the nullclines should be a straight line and a cubic. But all this "lumping" together is what I don't get exactly.

Problem 4.7.1

A few representative plots of what I got for 4.7.1. Not intended to be copied, just a reference.




Papers for Thursday 30th on Hydrophobic collapse in proteins:

  • Principles of protein folding -A perspective from simple exact models - Dill et al. - Protein Science (1995), 4:561-602
  • Initial hydrophobic collapse in the folding of barstar - Vishwas R. Agashe, M.C.R. Shastry, and Jayant B. Udgaonkar - Nature, v.377, 26 Oct. 1995

Papers for Thursday 30th on "continuous and discontinuous propagation in cardiac muscle":

The main paper for Tuesday, Nov 28 on the topic of the Phage Lambda can be found here: This paper's model is based on the model from Shea and Ackers (1985) which can be found here:

We can list our references and such here. Also, I found this site Max Anim which has some pretty good animations (available for download) for some of the processes we've seen in class, and that might apply to some of your projects (e.g. here is the Lac operon animation, which clarified a lot for me).

The first presentation is on Tuesday Nov. 21, on the subject of B12 riboswitches. The reference to look at is:

papers for our thursday talk