I do not understand how to implement the reaction mechanism to produce a ligand in a simple unimolecular reaction. I thought I'd want a two state reaction with equal rate constants and one state that produces a molecule of glutamate into the simulation. I can't figure it out from the MCell ref guide.

Here is MDL code that doesn't work

G_m = 5⁵ /* (/sec) transition that will make glutamate */

G_r = 5⁵ /* (/sec) return transition back to base state */

DEFINE_REACTION makeglu {

G[>GA{Gm:*glutamate, POSITIVE_POLE}]

GA[>G{Gr}]

REFERENCE STATE G {

glutamate NUMBER_BOUND = 0

}

}

I get probabilities for the transitions that seem reasonable, but I never actually change state in a run with a timestep of 1 microsec and # steps of 100000.

Can you clarify where I am messing up? I do not understand the need for the REFERENCE_STATE code, but found I got a parse error if I did not include it.

Your mechanism is almost there. What you need now is to add the stupid, undocumented "dummy ligand" step to your reaction (this step, and the REFERENCE_STATE will be unnecessary in MCell version 3, which will finally include molecule-molecule interations as well). But for now, here's how to use a dummy ligand to make a purely unimolecular reaction work properly:

DEFINE_LIGAND dummy {

DIFFUSION_CONSTANT = 0

}

DEFINE_REACTION makeglu {

R[>G{0:+dummy}]

G[>GA{Gm:*glutamate, POSITIVE_POLE}]

GA[>G{Gr}]

REFERENCE STATE G {

dummy NUMBER_BOUND = 0

glutamate NUMBER_BOUND = 0

}

}

What you do is define a make-believe ligand called "dummy" with any diffusion constant you like (diffusion constant doesn't matter because we'll never release "dummy" into the simulation. Then when you define your reaction you define a state that can bind "dummy" with any rate constant (again the rate constant doesn't matter because "dummy" will never be diffusing around to actually collide with the effector). Note that we define the binding of "dummy" to state makeglu.R which then transitions to state makeglu.G, so the state makeglu.G, by definition has one molecule of "dummy" stuck to it -- this is the essence of the trick. Now when you add effectors to your surfaces you add state makeglu.G --which will come with one molecule of "dummy" along for the ride. A subtlty of MCell's current reaction algorithm is that it's driven by the presence of bound ligand -- purely unimolecular reactions won't work correctly (this was a silly oversight in the original design of the algorithm). So, the "dummy" ligand trick is a work around for this oversight. Not a difficult work around but also not documented anywhere. We did go over it during the workshops but who's memory is that good?

You don't really need to have a transition to state GA if you don't want to. You could just do this:

G[>G{Gm:*glutamate, POSITIVE_POLE}]

This would allow each makeglu.G effector the chance to emit glutamate at rate Gm on every time step. This would simplify the math in figuring the rate of production of glutamate.

There will be a new way of iterating over molecules in the MCell-v3 with molecule-molecule interactions. Molecules which, in their current state, can diffuse or undergo any unimolecular state-change will be classified as "reactants" and will be checked for bimolecular or unimolecular transitions as appropriate. This scheme will require a slightly different MDL syntax but will no longer require the confusing REFERENCE_STATE statement. There will be numerous other advantages and our documentation is undergoing a MAJOR overhaul. All thanks to help and advocacy of colleagues like you, who have been instrumental in our success in landing new NIH funding for MCell development!

When I use dummy ligands, I don't even include the R>G transition. Instead, I just set the reference state number bound to one instead of zero for state G.

Can a reaction mechanism do the following:

a state change from A to B results in making 3 ligands of type x and 1 ligand of type y?

You can't do that all in one single reaction step but you can do it in several separate steps that transition to temporary intermediates. The new version of MCell with molecule-molecule interactions will allow what you're asking for in one step.

Here's one possibility to get you going with the current version of MCell:

A[>A1{K1:*x]

A1[>A2{K1:*x]

A2[>A3{K1:*x]

A3[>B{K2:*y]

As long as the rate of production of x and y is low compared to your time step this should work fine.

I've got the multicube arrangement similar to Franks, etal and I want to prevent molecules from going down the extracellular space paths between a small number of the cubes. I was thinking of creating a thin reflective box to cover the entry point into those paths, but I wanted to be sure not to violate the MCell rules on coincident surfaces. So, my question is, how accurate is MCell at separating surface locations? ie. what is the # of significant digits in coordinate specification? Would [0,0,0] versus [0,0,0.0001] be calculated as being different? If they are different, would molecules theoretically be able to make random walks in between these very closely apposed surfaces (albeit at very low probability)?

The only other way to achieve what I want would be to simply have a very large cube (and therefore no excell space) on the "presynaptic" side, but the way I am generating the array of cubes makes it impossible to do that. If the idea of the "almost coincident" box is not viable I will just brute force it.

On your MCell question. It really would help have to have picture of what you're thinking about when you say you "want to prevent molecules from going down the extracellular space paths between a small number of the cubes."

The rules on coincident surfaces are really guidelines to follow for building models. If you are very careful about what you're doing you can construct models with coincident surfaces in select locations to achieve a desired effect, such as blocking a diffusion path. I'll explain in more detail in a moment, but first let me say that the accuracy of surfaces in MCell is double precision and so separation between surfaces as small as 1 double precision epsilon (2.2204460492503131E-16) will look like a crack and particles will leak through it. Note that the units of this epsilon value is not in microns but is in MCell's internal length unit where one length unit is equal to 1/sqrt(EFFECTOR_GRID_DENSITY). If your EFFECTOR_GRID_DENSITY is 10000/um², one length unit will be 10nm so epsilon will be 2.2204460492503131E-24 meters (i.e. 1 billionth the size of a proton). This is pretty stinking small. To make it easier for MCell to detect and correctly handle coincident surfaces (if the model happens to contain any) MCell further limits the precision to 1e-12 length units (i.e. 1e-20 meters if the length unit is set to 10nm). Still pretty stinking small. The upshot of all this is that if you don't want particles to leak between surfaces they really have to be mathematically coincident. There are only two ways to guarantee that this is the case. 1) The polygons that make up a surface must share the appropriate vertices with one another. or 2) If the vertices are not shared then the polygons must lie in the same plane. Now with finite numerical precision, the only way that two polygons can lie in the same plane is if they lie in a plane parallel to one of the 3 cartesian coordinate planes (i.e. parallel to the XY, XZ, or YZ planes). For example, if two polygons have the same exact Z value for all of their vertices they lie in the same exact plane parallel to the XY plane. If you apply the same rotation transformation to these two polygons they will no longer lie in the same mathematically exact plane with each other due to rounding errors that arise from finite precision.

So to get to your situation, since your boxes are all created aligned with the cartesian coordinate system (guanranteed by the BOX statement as long as you don't ROTATE the boxes), you should be able to create proper coincident surfaces along, say, the TOP of some of your boxes as long as you create them at the exact same Z value as the TOP of the boxes. You'll also have to make sure that none of the coincident surfaces contain effectors and that all of these surfaces have the same surface permeability property (i.e. all REFLECTIVE, ABSORPTIVE, or TRANSPARENT). If you create your coincident surfaces this way things should work properly.

I have broken up my MDL code into logical/functional "modules" that I load as INCLUDE files after I've defined the time step, # iterations, ligand and diffusion rate.

I'd like to be able to modify a variable contained in one of the include files in the body of the main MDL file rather than having a group of different include files. I have tried this by having a line in the main MDL file that occurs after all the INCLUDE files have been specified, but that has not resulted in a parameter change.

Am I correct that I cannot modify the value of a variable set in an INCLUDE file (explicitly a rate in a kinetic model for an effector that makes ligand) by having a line of code in the main MDL?

You are allowed to modify the value of a variable at any time. The way you use the variable is the more important issue. Modifications to values of variables are not retro-active -- that is, if I do this:

• 1. a=1

• 2. b=a

• 3. a=2

• 4. c=a

• 5. a=3

printf("a = %g b = %g c = %g\n")

The output of the printf statement will be:

a = 3 b = 1 c = 2

So, if you assign a value to a variable in an INCLUDE file, then make use of that variable some how, then modify the value of that variable in your main file, and then finally make use of that variable again, your use of the variable will follow the behavior of the example code I gave above.

Does this make sense? Are you seeing something different from this behavior?

Just wanted to check which version of MCell is the correct one if I'm going to run on and SGI Altix with the Itanium (Madison) processors? I assume the Intel for Linux version will work since there are no other options, but I didn't know if it is an optimal version or not.

Yes, the linux x86 version is the optimal one. We have compiled MCell on the Itanium with gcc-3.2.3 and with Intel's optimized Itanium compiler with no significant change in benchmark results. In addition, the x86 binary compiled on an x86 machine runs as fast on the Itanium as the binary built on the Itanium. Bottom line: the Itanium is not a worthwhile MCell platform. Now, we've also tried MCell on the AMD Opteron and wow! The 1.5GHz Opteron we tried was 10% faster than a 3GHz P4 Xeon on our MCell benchmark.

The MN supercomputer institute folk are trying to get me time on a variety of their resources. For the next couple of weeks their SP is still running AIX 4. I think the binary for download from PSC still has the IP lock in the code as I got the MCell Licensing Error when I tried to run using that version. I haven't had that problem with the AIX 5 code. I assume the AIX 4 is considered somewhat obsolete and so, has not been recompiled. Is that a lot of work?

If I can get the no IP lock version, I can benchmark my runtimes and then use that to lobby for machine time that might make more sense from the standpoint of turnaround times, etc.

Just following up on the SGI Altix. I'm seeing this machine being twice as slow as the Linux cluster which has 1.2 GHz PIII processors. I did not expect the performance to be so lousy and I wanted to check if you thought it was possible. As far as I can tell, I ran the identical MDL files on both systems. This is relevant, as I can get lots of processors on the 2 Altix systems, but it may not be worth a lot if the runtimes are, comparatively, abysmal.

The AIX 4 code has not been upgraded because all the SP machines we have access to have been upgraded to AIX 5.1 so we no longer have access to any AIX 4.X machines, as far as I know. Rebuilding MCell only takes a few minutes so it's not a big deal. If I had an account on one of your AIX 4 machines I could do it very quickly.

On the Altix systems I'm assuming it's a 1.5GHz Itanium2 machine, correct? If so I'm not surprised by the benchmark because the binary you're using is a IA32 binary not an IA64 binary. You can verify this by typing:

file mcell

This will report all sorts of goodies about the mcell binary you're using. Now, the Itanium2 (which is an IA64 platform) can run IA32 binaries but in a highly degraded performance mode. To get full performance you'll need an MCell binary built for the IA64. Then you'll see performance on the 1.5GHz Itanium2 which is about 90% the speed of a 3GHz P4 Xeon. If you had an AMD Opteron machine (which is an AMD64 platform) \ the performance of binaries built for IA32 and IA64 would run equally fast and, for a 1.4GHz Opteron, would be about 110% the speed of a 3GHz P4 Xeon. This is why Los Alamos and Lawrence Livermore have gone with the Opteron for their new clusters. (Yeah, great, now they tell me) :-)

In any case, If I also had an account on your Altix I could build MCell for IA64 there, as well.

I see that OpenDx is up to v 4.3. I'm still running 4.0. Are there good reasons to upgrade or is it fine just to stay with what I have.

Rule #1 of software upgrades: If it ain't broke don't fix it.

I have gotten the following statement from MCell being run on an IBM Regatta.

"MCell: error on line: xxx of file: xxx.mdl cannot store subvolume partition data."

The only difference between a run that had no error and this run was that the world was populated by smaller cubes (.25 micron side length) when the error occurs versus 0.5 micron side length when it does not. The partitioning statements are identical.

I had specified 4 Gig of RAM for this job.

Is this just a case of needing more RAM? I won't be able to get another job through the queue for a day or more and I thought I'd check with you in case it might be something else.

The same MDL file that gave the error on the Regatta did not give the error when run on an interactive node on a PIII linux cluster.

Yes, it's an out of memory error message. Any message from MCell that says it "cannot store" some item means you've run out of memory. Keep in mind that the number of required partitions (and associated memory usage) scales with the cube of length of sides of the partitions. By going from 0.05 microns (I assume you meant to say 0.05 micron partitions) to 0.025 microns the number of partitions and memory usage would have increased by a factor of 8.

I believe the IBM Regatta is a 64bit architecture and so uses 64bit pointers -- as compared to 32bit pointers on 32bit architecture machines such as the Intel Pentium. This will make the IBM Regatta consume much more memory to compute the same MCell simulation since MCell makes extensive use of data structures that contain pointers. This is to be expected and is part of the cost of 64bit computing. An analogous price was paid in the migration from 16bit to 32bit computing.

I do have an mcell checkpoint question. Am I correct that I cannot change the diffusion rate of molecules between checkpoint runs?

You can change the diffusion rate of molecules between checkpoint runs.