With increased computing power, computer simulations may be performed for larger-scale models. Computer simulation may be used to help solving problems in a very efficient and economic way, and testing or verifying scientific hypotheses so as to better direct research experiments. All simulations essentially encompass three steps: a. modeling; b. simulation; c. visualization (or result analysis). These three steps are strongly correlated: While a well designed model is key to success, visualization also plays an important role: since simulation data can be converted into 3D graphics, it is possible to compare a model with the real object and as a result further improve the simulation. Table 1.1 shows the basic steps involved in a computer simulation.

Table 1.1 A Generalized Computer Simulation

Modeling is the first step for any kind of simulation. A model is the definition of an application and a detail plan about how to perform this application. A well designed simulation model should include all the information needed by the simulator, such as simulation objects (physical objects and logical objects), simulation parameters, simulation mechanisms, optimizing algorithms and so on. The model may be translated into a simulation input file in the format required by the simulator to perform the simulation.

MCell was initially designed for microphysiological simulations. To simulate a realistic microphysiological system, representations of the relevant cellular and sub-cellular structures must be generated at very high resolutions (the electron microscope level). To achieve this, we first use biochemical methods to treat some tissue for Electron Microscope (EM) to take pictures (this step is also called image processing). Then the EM image is scanned to transform the image into digital data, the digital data are then used to reconstruct a 3D digital image for further study (simulation); this step is called 3D reconstruction. The simulation object may be an object from the real world ( part of neural muscular junction, a channel, a clamp, a piece of membrane etc.) or a computer-generated three-dimensional object. The user may also use the MCell Model Descriptive Language (MDL) to define a physical object as the 3D simulation object.

After you have the simulation object and the simulation model in hand, you are ready to run a simulation. The simulation will be done in three steps:

• Simulation initialization;
• Time-step iteration;
• Output the simulation results.

After a simulation is finished, the user may collect the data for visualization; the simulation model may be further improved based on the effect of the visualization. A more refined simulation can be started after modifications of the model.

MCell was initially designed for computational neurobiology but it is not limited to that area. It can also be applied to a wide range of micro-studies such as diffusion-reaction applications, molecular dynamics, and computational chemistry. The MCell team is working to add electronic physiology into the simulation in future versions.

In computational neurobiology, there are four typical events that are modeled by an MCell simulation. They are:

• The release of molecules from a structure (e.g., vesicle release);
• Molecular diffusion within spaces defined by arbitrary surfaces (e.g., pre- and postsynaptic membranes, or a cell membrane with attached patch clamp micropipette);
• Creation or destruction of molecules (e.g., synthesis, hydrolysis processes, or redox reactions);
• Chemical reactions undergone by ligand and "effector" (e.g., receptor or enzyme) molecules.

Ligands, effectors, reaction mechanisms, 3D surfaces, and other simulation components are specified using a simple descriptive language in MCell. The language is called MCell Model Description Language (MDL) and was designed with the biology-oriented users in mind - it is used to describe the simulation model.

MCell can be executed under different operating systems: Windows, Windows NT, Free BSD, and UNIX. The MCell execution command is the same for each operating system. If the simulation input file is file.mdl, then the MCell execution command will be:

mcell [options] file.mdl

The MCell simulation input file is defined by the MCell Model Description Language (MDL), and it provides MCell with all the necessary information for a simulation, whcih may include simulation components, simulation objects, output modes, simulation mechanisms and the simulation structure. There are no restrictions on the name of the input file, the file extension .mdl is not required by the MCell. We recommend the name format file.mdl only for the convenience of file organization.

First, one or more MDL input files are interpreted to create the simulation objects. Next, MCell will start the simulation iterations. Each MCell iteration equals one Monte Carlo time-step. Via a feature called checkpointing, MCell can break up a simulation into smaller sets, each one simulating the model over given amounts of time. The simulation is based on the following assumptions being true at each time step:

• Every ligand (small molecule) has a position and is flagged as either free or bound by an identity value (flag).
• Free ligands move randomly and may interact with effectors (large molecule).
• The position of each effector is fixed; each effector is considered to be in one of the states of a given reaction schema.
• The state of an effectors is associated with a surface area and a probability of binding, given that a transmitter molecule hits the effector surface.

MDL files are order-dependent. Whenever your run MCell with an MDL input file, MCell will interpret the MDL file in a top-to-bottom order. Initialization will be done during MDL parsing. Most MDL keyword statements are optional, which means you only use them as needed. However, there are two required MDL keywords for all MDL files: time step of each simulation iteration and the total number of iterations of the simulation. The anatomy of a well organized MDL file is given below.

Figure 2.1 General anatomy of an MDL file. MDL files are order dependent as shown in the figure. Comments may appear anywhere inside these files; time step and number of iterations are the only required parameters.

You may assign text variables,numerical variables, and numerical arrays (whose elements are numerical values or numerical expressions) in a MDL file. You may also express text and numerical expressions in MCell. A text expression is a series of concatenated text, and a numerical expression is an arithmetic relational expression.

MDL recognizes text variables when they are enclosed by quotation marks. All numerical data are regarded as double precision internally, so there is no need for variable-type definition when creating numerical variables. Also, user-created variables do not need to be capitalized like MDL keywords, and they may be used anywhere in a MDL file. There is no restriction for naming variables, but by using carefully chosen variable names you can make your program self-documenting - MDL lines themselves would tell whoever might be reading the code what it does. By reading the defintion in Table 2.2-1, you will learn how to create your own simulation variable.

Table 2.2-1 How To Create The User Variables

MDL supports the typical arithmetic operators and numerical operations. Arithmetic operators can apply to both the user-created variables and MDL-defined variables (the MDL keywords). Numerical expressions can appear anywhere inside or between MDL blocks in a MDL file.

MCell supports five standard numerical operators, they are:

A numerical expression combines standard arithmetic operators (+, -, *, /, ^), mathematical functions, the value p (can be expressed as PI in the MCell), user created variables, and MDL keywords. The available mathematical functions are given in Table 2.3-1.

Table 2.3-1 Mathematical Functions of MCell

Like any other programming language, comments are widely used in MDL files to increase the readability of the files. Comments are allowed anywhere (inside or outside the functional block) in a MDL file. The MCell parser will ignore comments during simulation.

Comments are any contents stated between the start tag (/*) and end tag (*/). The format of a MCell comments is shown in Table 3.1-1.

Table 3.1-1 The MCell Comments

In an MCell simulation, molecular diffusion and chemical reactions are the two major applications. Parameter time plays an important role in both by controlling the simulation process. Total time of a simulation is the product of number of iterations and time step of the iteration. The MDL commands ITERATIONS(total number of iteration of the simulation) and TIME_STEP(time interval between two adjacent iterations) are the only two required simulation parameters in a MCell simulation. The simulation will print out an error message if you try to run a simulation without specifying ITERATIONS and TIME_STEP.

Table 3.2-1 shows the definition of the MDL command ITERATIONS and TIME_STEP.

Table 3.2-1 ITERATIONS and TIME_STEP

The MDL command ITERATIONS defines the number of iterations that will be executed in a MCell simulation. There is no default value for this parameter in the MCell simulator. It is a required simulation parameter for any MCell simulation.

The MDL command TIME_STEP specifies the time interval between each iteration. TIME_STEP is expressed in seconds and there is no default value in the MCell simulator for this parameter. TIME_STEP is a required simulation parameter. TIME_STEP is also the only parameter that would effect the convergence of a MCell simulation.

INCLUDE_FILE is the MDL command to include an existing MDL file into a simulation input file. It can only include one file per usage, but you are allowed to use the command INCLUDE_FILE as many time as you like in your MDL file. There are two ways to include files using the command INCLUDE_FILE:

• Double quote the name of the file that you are going to include and assign it to the keyword INCLUDE_FILE;

  INCLUDE_FILE = "file_name"

• If the included file name has been assigned to a user-defined variable in the pre-included file, you can assign the related user-defined variable name to INCLUDE_FILE.

  INCLUDE_FILE = variable_name

There are no naming restrictions for the included file. When you assign a filename to INCLUDE_FILE, just double quote the file name; No double quotes are needed for a variable name. You can freely use variables and parameters defined in the included file after the INCLUDE_FILE statement in the simulation. Figure 3.1 shows the usage of the command INCLUDE_FILE in an MCell simulation.

Figure 3.1 If you have the function blocks and simulation parameters saved in two separate MDL files (as File1.mdl, File2.mdl in the figure), you can use the MDL command INCLUDE_FILE to add them into your simulation input file. You can have as many include files as you need.

We give two examples to introduce how to use the MDL command INCLUDE_FILE in a MCell simulation.

It is possible to insert check points to check or adjust a simulation before it terminates. This method is called checkpointing. Checkpointing is a commonly-used method to control a long simulation process; it is also a powerful and simple way to change run-time simulation parameters.

The checkpointing method is performed by inserting checkpoints into specific iterations of a simulation. When the simulation reaches a checkpoint, it will pause at that point. Checkpoints, therefore, can split a long simulation into several segments that can be executed sequentially. If you have M checkpoints, there will be M+1 simulation segments. You may use a Perl script file to operate the simulation automatically, the script file can load each simulation segment sequentially after each checkpoint. The three MDL keywords used to define the MCell checkpoints are:

• CHECKPOINT_ITERATIONS
• CHECKPOINT_INFILE
• CHECKPOINT_OUTFILE

Notice: The MDL command CHECKPOINT_ITERATIONS is a required argument for the checkpointing method, the other two are optional.

When a simulation pauses at a checkpoint you can resume the simulation from that checkpoint by assigning the name of the CHECKPOINT_OUTFILE (which is defined in the previous simulation segment) to the MDL command CHECKPOINT_INFILE in the following segment. The necessary condition to continue the MCell simulation at a checkpoint is that the CHECKPOINT_OUTFILE is defined at the checkpoint, so that the simulator would store the previous segment to the checkpoint, and later use that information to resume the simulation. CHECKPOINT_OUTFILE will be read into the next segment by CHECKPOINT_INFILE. If CHECKPOINT_INFILE is not defined in the following segment, the simulation will start over from the beginning. Definition of checkpointing is given in Table 3.4-1.

Table 3.4-1 Definition of The Checkpoint

The MDL command CHECKPOINT_ITERATIONS specifies the location of the checkpoint within a simulation. The MCell simulation will automatically stop when it reaches the point of CHECKPOINT_ITERATIONS. The checkpoint method is only valid when CHECKPOINT_ITERATIONS is less than ITERATIONS.

The CHECKPOINT_OUTFILE saves the simulation results before the checkpoint. CHECKPOINT_OUTFILE should be defined after CHECKPOINT_ITERATIONS. After the checkpoint is reached, the primary simulation results will be saved to the file specified by CHECKPOINT_OUTFILE. These results may be read by the command CHECKPOINT_INFILE in the next simulation segment.

The simulation input file is split into sequential segments by checkpoints. In order to maintain consistency, the MDL command CHECKPOINT_INFILE should be specified in the second simulation segment. The output file name defined by CHECKPOINT_OUTFILE in the previous segment should be assigned to CHECKPOINT_INFILE in the next simulation input segment. The command CHECKPOINT_INFILE is an optional statement for the checkpointing method; however, if no CHECKPOINT_INFILE is specified, the simulation will begin from the beginning of the simulation, iteration 0.

A MCell simulation is a general Monte Carlo simulation. Like the Monte Carlo simulator, MCell focuses on molecular diffusions of short distances. The majority of movements of the simulated molecules is diffusion, which is an irregular movement that will result in a uniform distribution of the diffusing molecule from the initial non-uniform state. This type of movement is characterized as Brownian motion; therefore the movements of the molecules are simulated by an optimized Brownian algorithm

Because of the random characteristic of molecular diffusion, the step length and the moving direction of a molecule at each iteration time step are chosen randomly. The diffusion equation is solved first to find the key parameters of molecular diffusion. For computational conveniences, movements of molecules are expressed in polar coordinates. If we denote the concentration of the simulated molecules C(r, t) as a function of the radial distance r, and time t, with a diffusion constant D_l, then the diffusion equation expressed in polar coordinates will be an one dimensional equation:

Equation (3.1) may be solved analytically. If the total amount of molecules in the whole system is M, then the concentration of the molecule C(r, t) will be expressed as:

Since the total number of molecules (expressed as N) in volume dV is the product of the volume and the concentration of the molecule in the system. Then N will be computed by:

We define p_r as the fractional probability. It is the probability of obtaining a net radial displacement between distance r and r+dr for each single diffusing molecule. The volume between distance r and r+dr is dV, number of molecules in dV is N, so the the fractional probability p_r is N divided by the total number of molecules in the whole system, which is:

We use the fractional probability p_r to weigh the movements of molecules. After weighing, the mean radial displacement (l_r: the step length of the molecule diffusion.) of the molecule would be the expected value for a radial distance r. If we integrate p_r over r for each molecule at each time step to get the mean radial step length l_r, the computation will be too expensive. The MCell simulator optimizes this computation by introducing a dimensionless parameter s into the computation, s is defined as:

Substitute s into equation (3.4), we get:

By integrating equation (3.6), we obtain a set of s with equal probability values. If we substitute s into equation (3.5), we will get a set of radial step lenghts with equal probability values for any molecular diffusion. The radial step length is the product of s and the constant (4D_ldt)^1/2. You will have as many radial step lengths to choose from as the number of s values you have. Since all radial step lengths have the same probability, and diffusion of the molecules are random, we only need to pick one radial step length from the set randomly to move the molecules at each time step. MCell's optimization is to save all the equal probability s values into a look-up table. During each MCell time step iteration, we only need to pick one s value randomly from the table, multiply s by the constant (4D_ldt)^1/2 to get the radial step length. This MCell optimized algorithm is faster and easier compared to the traditional Brownian algorithm, which integrates the probability value for every time step to compute the radial step length.

The more s values we have, the more choices we have for the radial step length. Precision of the simulation therefore increases with the number of s values available. The trade-off is that more computer memory is needed to save the look-up table for s. You may assign the number of s by the MDL command RADIAL_SUBDIVISIONS. The default value of the RADIAL_SUBDIVISIONS is 1024 and the maximum value of s is 4096.

Knowing the radial step length and radial direction, molecules are moved from one point P1 to the next point P2 using ray tracing algorithms. Ray tracing was originally a technique used for rendering realistic images of 3D scenes in computer graphics. Imaginary, where rays are extended against the objects in the scene and the closest intersecting object is found.

In each MCell time step, we imagine that a ray is extended from the original location of each molecule to the selected radial direction. The ray will go through a distance determined by the computed radial step length, then we determine the location of the point P2 and any intersection(s) of the ray with the objects on its way.

If the ray doesn't hit anything, it means that the molecule will not intersect with any object as a result of its movement, so we move the molecule from point P1 to point P2 along the picked radial direction with the computed step length.

In the case that the ray does intersect with an object, the first step is to judge what kind of object it is; we should know if it is a wall of a physical object or an other molecule. If it is a wall of a physical object, we need to check further if there is a molecule on the specific surface area. If there is no molecule, we then continue to check the surface permeability of the wall:

• A reflective wall would reflect the ray and a new destination would be computed;
• If the wall is transparent, it would let the molecule go through and finish its trip to the destination P2;
• The wall could also be absorptive, in which case the molecule will be absorbed by the wall, removing it from the molecule list.

RADIAL_SUBDIVISIONS is the number of divisions in an unit length along a radial direction. Therefore the RADIAL_SUBDIVISIONS are a set of step lengths with equal probability values. These equal probability values are saved in a lookup table during MCell initialization; they are used to compute the radial step length of the molecule during the time step iterations. The default number of RADIAL_SUBDIVISIONS is 1024, which provides enough precision for the general MCell simulations; larger RADIAL_SUBDIVISIONS value has a trade off of larger computer memory requirement.

Table 3.5-1 The MDL Command RADIAL_SUBDIVIONS

RADIAL_DIRECTIONS is the number of the radial directions in a single octant of a unit sphere. They are the possible radial directions that we can use as the molecules' diffusing directions. The directions are saved into the second lookup table during MCell initialization. The default RADIAL_DIRECTIONS is 16384, which requires approximately 400 Kbytes computer memory for the computation, and the maximum value will need almost 3.2 Mbytes memory.

Table 3.5-2 The MDL Command RADIAL_DIRECTIONS

EFFECTOR_GRID_DENSITY is the number of the effector grids per unit area on the surface of a simulation object. The effector grid is the mesh element on which the effector sites could be placed. This value decides the number of the effector tiles that will be distributed on the surface of the simulation objects. The effector tiles are possible locations for the effector sites. The default EFFECTOR_GRID_DENSITY value is 1 mm^-2. The definition of the MDL command EFFECTOR_GRID_DENSITY is given in Table 3.5-3.

Table 3.5-3 The MDL Command EFFECTOR_GRID_DENSITY

Spatial partitioning is the MCell run time optimization. Spatial partitions are transparent planes that the user places on the simulation object along the X, Y or Z axes to make spatial sub-volumes. To further understand the advantages of spatial partitions in an MCell simulation, we introduce the computational mechanism of the MCell simulator. All MCell simulations are carried on spatial objects which occupy certain volumes in space. From the computational point of view, the movements of the molecules inside the simulation objects should be calculated with the consideration of the whole volume, we call this volume the computational volume. If we have a larger number of molecules moving inside a bigger object, it means that we need to simulate more molecules inside a bigger computational volume at each time step and it will require more time and more memory. The larger the object is, the more computer memory needed. If we use spatial partitions to divide the original spatial volume into smaller sub-volumes, then the size of computational volume is reduced to that of each sub-volume; at the same time the number of the molecules inside each computational volume is reduced to those inside each sub-volume. As a result, the size of the computational object and the number of molecules decreases substantially; consequently the computation for the intersections between the molecules, as well as the computation for the movements of the molecules are reduced. Thus a well designed spatial partition will increase the simulation speed dramatically.

Table 3.5-4 The Spatial Partitioning in The MCell

There are two types of defined objects in a MCell simulation: logical objects and physical objects. Logical objects are objects without definite or fixed physical locations in 3D space. Objects like ligands, ligand release patterns and chemical reaction mechanisms are logical objects defined in the MCell simulator. Physical objects are objects with definite locations in 3D space. In this section, we discuss the definition of logical objects in a MCell simulation.

A MDL function consists of a function name (capitalized MDL keyword) and a function body (contained within the brace {}). Each MDL function contains both required and optional elements.

Table 3.6-1 Definition of The Molecule

The diffusion constant of the defined ligand. This parameter measures the rate at which particles can diffuse through a medium; it depends on the substance, the medium, and the temperature. The diffusion constant has a dimension of area/time. In a MCell simulation, this value is used to calculate the radial step length of the ligand.

The reference diffusion constant. This value is used to calculate the probability of ligand binding.

Charge is the electrical charge that a ligand carries. The movements of the ligand form a flow of electrical charges, called current flow. The parameter CHARGE is particularly useful for the simulation of ion channel.

Sometimes you may want to define a temporal pattern to control the ligand release process. The release pattern is a series of discrete events, all of which are defined using the simulation start time t = 0 as reference. Each release is a pulse as illustrated in Figure 3.2.

Figure 3.2 Timing Pattern of the ligand releasing

The temporal pattern of ligand release shown in Figure 3.2 is typically defined as "a train". The pattern of train distribution is the pattern of the ligand release. The release time depends on the number of trains, train duration, train interval and delay time. Each train contains one or more release events (pulses). The release pattern depends on the train duration, release interval and the train interval. The ligand release pattern is defined by the MDL function DEFINE_RELEASE_PATTERN with a user-given name. The defined release pattern will be assigned to the MDL command RELEASE_PATTERN in the function SPHERICAL_RELEASE_SITE by the defined pattern name. Table 3.6-2 shows the format of the function DEFINE_RELEASE_PATTERN.

Table 3.6-2 The MDL Function DEFINE_RELEASE_PATTERN

The time interval between the simulation start time t=0 and the first ligand release event in the ligand release process.

The time interval between two consecutive release events (pulses) within a release train. It is the inverse of the frequency of the release event.

The time interval between the beginning of one train to the beginning of next. It is the inverse of the train frequency.

The duration of each train. It is the product of the number of events in a train and the RELEASE_INTERVAL of the train.

Table 3.6-3 Chemical Reactions in MCell

We give two examples of the DEFINE_REACTION function: the single path reaction and the multi-path reaction, shown in Figures 3.3 and 3.4, respectively. In Figure 3.3, we give a single path chemical reaction, in which state A will react with (bind) state E when the reaction constant k1 is reached, this reaction will produce a new state--AE; State E and A will unbind and be released from state AE when the reaction constant k2 is satisfied.

Figure 3.3 Single path chemical reaction

We can write this reaction with the MDL function DEFINE_REACTION with two simple statements as in the following example.

We express a multiple path chemical reaction path by path, and each path is enclosed in a separate set of square brackets. For the reaction shown in the Figure 3.4, we have four paths:

• First, state A binds with state R⁰ as the reaction constant k11 is achieved; this reaction will produce the state AR¹;
• Second, state AR¹ will unbind and release both states R⁰ and A as the reaction constant k12 is satisfied;
• Third, state A binds with state R⁰ when the reaction constant k22 is reached, and state AR² will be the new state that is formed after the reaction;
• Fourth, state AR² unbind and release both states R⁰ and A as the reaction constant k21 is satisfied.

Figure 3.4 The Multiple path chemical reaction

The MDL expression of this reaction is given in the following example:

The ligand interaction operators are the symbols used to tell the MCell simulator what kind of interaction will happen between the defined effector state and the selected ligand molecule. The symbol % that appeared in the function DEFINE_REACTION represents these 6 ligand operators (six kinds of the ligand interactions):

Operator "+" is the ligand binding operator; In the function DEFINE_REACTION, when this operator appears, the ligand behind it will bind with the defined effector, under the proper polarity conditions.

Operator "-" is the ligand unbinding operator; In the function DEFINE_REACTION, when this operator appears, the ligand behind it will unbind with the defined effector, under the proper polarity conditions.

Operator "~" is the ligand transport operator. The ligand which follows this operator will be transported through the medium. An example would be molecular movements through a transmemebrane channel.

Operator "#" is the ligand destruction operator; In the function DEFINE_REACTION, when this operator appears, the ligand behind it will be degraded by the defined effector when the chemical rate constant is satisfied.

Operator "*" is the operator for the production of the diffusing ligand.

Operator "@" is another ligand producing operator, the ligand production of this operator has a Poisson distribution, so we call it Producing Poisson.

Besides ligand interaction operators, you also need the ligand polarity values to define the chemical reaction. The MDL keyword polarity is used to define the binding and unbinding between the ligands and the effectors in the function DEFINE_REACTION. The value of polarity depends on the ligand interaction operators. Table 3.6-4 lists the polarity values for different ligand interaction operators.

Table 3.6-4 The polarity of A Ligand

Physical objects are objects with definite physical locations in 3D space. Physical objects may be objects created from 3D reconstruction, objects designed by computer graphics tools, or objects defined by the MCell physical object definition. They are designed to hold or locate the logical objects in MCell simulations. All simulation events occur either on the surface of physical objects or around them. Three basic physical objects defined in the MCell are:

• The ligand release site ;
• The box object;
• The polygon list object;

A ligand release site is the physical location of a ligand source in 3D space. The MDL function SPHERICAL_RELEASE_SITE is used to define a ligand release site. This function provides the MCell simulator information about the location and size of the site, the number of the ligands that are released, the name of the released ligand, and the temporal pattern of the ligand release process. The ligand release source is modelled as a sphere by the MDL function SPHERICAL_RELEASE_SITE. You can create a point source site by assigning zero to the diameter of the ligand source.

The reason for this simplification is that the ligand source is usually very small in the physical world, and most of them are round shaped. So for computational convenience, we neglect the shape of the ligand site and regard it as a small sphere. But a sphere is not the only shape that may be used for the ligand source, you may define an arbitrarily shaped ligand source too. To define an arbitrarily shaped ligand source, you first define a point source using the MDL function SPHERICAL_RELEASE_SITE and put it inside the desired shaped physical object. This physical object will be the ligand source. The arbitrarily shaped physical object may be defined by the MDL physical object definitions and other 3D reconstruction and graphical tools.

The MDL function SPHERICAL_RELEASE_SITE defines the ligand release site. Since we regard it as a spherical source, we only need to provide the origin and the diameter of the sphere in the function SPHERICAL_RELEASE_SITE to locate this source. The ligand source site should to have an uniform distribution of the ligand molecules. Geometric transformations, which include object translation, rotation, and object scaling, can be applied to the spherical release site as well as all the other physical objects. The Format of the function SPHERICAL_RELEASE_SITE is shown in Table 3.7-1.

Table 3.7-1 The MDL Function SPHERICAL_RELEASE_SITE

The x-y-z coordinates of the sphere origin defined in the Cartesian Coordinate system.

The total number of the selected ligand molecules that will be released from this source. There are three different ways to define the ligand release number.

The probability of the ligand being released from the ligand source. The default value of this parameter is 1, which means that the ligand will be released by default from the defined ligand release site. You may use this value to limit the number of ligand to be released.

The diameter of the spherical ligand source.

We have mentioned that the user can define a ligand release site with an arbitrary shape. You can follow the following steps to define an arbitrarily shaped ligand release site in your simulation.

1. 1.First, construct a physical object with the desired shape. The physical object can be defined from 3D reconstruction, built by graphical tools, or built by the MDL physical object definitions.
2. 2.The second step is to set a point ligand source using the function SPHERICAL_RELEASE_SITE, You can do this by setting the keyword SITE_DIAMETER in the function SPHERICAL_RELEASE_SITE equal to zero. The most important step is to locate the point source inside the designed ligand release site--the physical object that you defined in the first step.
3. 3.The third step is to define all the surface elements of the ligand release site as reflective, so that the ligand can diffuse and mix well inside. After a certain number of iterations the distribution of ligands inside the site will be uniform.
4. 4.Set a checkpoint in the simulation such that when the checkpoint is reached, the ligands inside the release site would have been distributed homogeneously. After the checkpoint, change the permeability of the release site surface elements to transparent, so when the simulation start again the ligands can escape from your designed ligand source.

The box object is a hexahedron constructed with the constraint that each face of the hexahedron is perpendicular to one of the x-y-z axes. We define the box object by specifying two vital points: the lower-left-front (LLF) corner and the upper-right-back (URB) corner (as point P1: lower-left-front point (LLF) and P2: upper-right-back (URB) point in the Figure 3.5). The box object is defined by the right-hand rule in the Cartesian coordinate system. To define the physical objects, we indicate whether the object is a closed object or an opened object.

Figure 3.5 The Box object defined in the MCell Simulation. The coordinates of the points LLF point and URB point should be given by user to define the box object in the MCell simulation.

Table 3.7-2 The MCell Box Object

The lower-left-front point (p1 in the Figure 3.5) and the upper-right-back point (P2 in the Figure 3.5) are the corner points of a box object which we use to define the location of it in a MCell simulation. The coordinates of the points P1 and P2 are the values assigned to the MDL keyword CORNERS. The x-y-z coordinates of the points P1 and P2 should be contained within the square bracket [], and separated by a comma.

The MDL keyword FULLY_CLOSED is used to describe the closed or open status of the defined object. Boolean values YES, NO, TRUE,and FALSE will be assigned to the keyword FULLY_CLOSED. The boolean value YES and TRUE have the same meaning, which is that the defined object is a closed object, boolean value NO and FALSE mean that the defined object is an open object.

The x-y-z vertex coordinates that we will use in the definition of the polygon list object are defined in the right-hand Cartesian coordinates system. The sequential numbers appeared in the function VERTEX_LIST are the index number of each vertex, and we will use these numbers to define the connections between each vertex by the function ELEMENT_CONNECTIONS. The definition and format of the MDL function POLYGON_LIST are shown in Table 3.7-3.

Table 3.7-3 The MDL Function POLYGON_LIST

The method to define a meta object in the MCell is given in Table 3.7-4. The function OBJECT has multiple applications in the MDL:

• The function OBJECT may be used to define meta objects;
• The function OBJECT may be used to copy a pre-defined physical object into meta object;
• The function OBJECT may work together with the MDL keyword INSTANTIATE to instantiate physical objects as simulation components.

Table 3.7-4 Definition of The Meta Object

Table 3.7-5 Creating The Simulation Components

The name of the instantiated object will be used as the domain name in the output specification. The simulation results of the physical objects can be specified by its domain name and its own object name which is defined in the instantiated object.

In MCell, you are able to define a desired surface region on an object template, and place effectors on the defined regions.

The term surface_modifier applies to all operations applicable to the surface elements of physical objects. surface_modifier includes defining and removing the surface elements, changing permeability of the surface elements, and adding effectors to the surface of the physical objects. You can specify an exact number of effectors or you can specify density of effectors on a given surface region.

The surface elements of the MCell physical objects are defined by the MDL command ELEMENT. Each value that is assigned to the keyword ELEMENT is called an element_specifier. The element_specifier varies with different physical objects. Besides defining the surface elements, you may also remove some surface elements by the MDL command REMOVE_ELEMENT. The specified surface elements are assigned to REMOVE_ELEMENT by the element_specifier. One ELEMENT, or REMOVE_ELEMENT statement can only specify one surface element. The format of the element operations are:

ELEMENT = element_specifier

REMOVE_ELEMENT = element_specifier

Different cases of the element_specifier are listed in Table 3.8-1.

Table 3.8-1 ELEMENT and element_specifier

From the Table 3.8-1 you may find that the command ELEMENT has different selection when it is applied to the box object and the polygon list object:

Permeability is an optional sub-function in physical object definitions. In order to define permeability, you should provide the name of the permeable ligand and the specific surface elements of the physical object. The format of the permeability definition is:

If a surface element is defined as transparent (MDL keyword: TRANSPARENT) to some specified molecule, then that molecule can go through the surface element as if there is nothing there.

An absorptive (MDL keyword ABSORPTIVE) surface will absorb the specified molecule. When the specified molecule hits the absorptive surface, it will disappear from the simulation forever.

If a surface element is reflective (MDL keyword REFLECTIVE), the molecule's direction will be changed after it hits the surface.

A new feature of MCell lets the user define the surface region on object templates; the user may place effectors on specific regions of the simulated object. This way the user can control the distribution of the effectors based on a realistic model.

Two general applications of this function are:

• Making a template object, defining a region on that object, and placing effectors on that region;
• Making two copies of the template, adding a second region to one copy, placing effectors on both regions;

Table 3.8-2 The MDL Function DEFINE_SURFACE_REGIONS

The effectors are large molecules used in MCell simulations, such as receptor proteins on the plasma membrane, transporter proteins, enzymes on an intra- or extra-cellular scaffold etc., it is also called the effector site. The effector site can bind and unbind with the surface of the physical objects.

There are three different ways to add effectors to the physical object:

• Adding effectors to the surface of the whole simulation object: the MCell simulator will distribute the effector over the surface randomly according to the EFFECTOR_GRID_DENSITY and the DENSITY of the effectors;
• Placing an exact amount effectors over the defined surface regions;
• Placing the effectors over the defined surface regions based on the specified effector density.

Table 3.8-3 The MDL Function ADD_EFFECTOR

The state name of the added effector. The value of STATE is one of the state_name defined in the chemical reaction specification.

The number of the effector molecules over an unit area. The effector molecules will be distributed on the effector tiles according to this value. The effector tiles are the triangular grids which divides the mesh elements on the physical object surface. The number of the effector tiles created depends on the area of the mesh element and the value of EFFECTOR_GRID_DENSITY. MCell divides the mesh element of the physical object into the effector tiles using the barycentric subdivision method. Each effector is an effector tile that has a chemical reaction mechanism associated with it, the effector sites are chosen from the effector tiles by a random number.

The MDL command ELEMENT is the surface element specifier of the physical objects. The command ELEMENT has different values for different types of physical objects. All the surface elements assigned to ELEMENT are represented by the symbol element_specifier in this menu. A detailed description of the command ELEMENT and its application will be introduced later in this section .

The command POLE_ORIENTATION is the orientation of the added effector sites relative to the normal vector of the physical object surface. Figure 3.6 illustrates the normal vector of the physical object surface. The face of the object that the normal vector leaves is the front face of the surface mesh element. The face that the normal vector enters the back surface. The POLE_ORIENTATION specifies which side of the object surface that the added effector lands on.

MCell provides new application for the distribution of the effectors. It can place effectors on defined surafce regions, either by an exact number of effectors or by effector density. To add effectors to the defined surface regions with an exact number, please follow the defintion given in Table 3.8-4.

Table 3.8-4 Adding an exact number effectors to the defined regions

To add effectors on the defined surface region by density is similar to what we did to add the exact number of effectors. What we need to do is to replace the NUMBER statement with the DESNITY statement. It is defined in Table 3.8-5.

Table 3.8-5 Adding effectors to the defined regions by density

The MDL function TRANSLATE is the operation used to move the physical object by a distance d_x along the X axis, a distance d_y along the Y axis, and a distance d_z along the Z axis. The object is moved by translating each point of the object to fulfill this operation. If point P(x, y, z) is a point on the physical object, after translation it will become to the new point P'(x', y', z'). The coordinates of the new point after the movement can be calculated from:

For computational convenience, we express the points P and P' as column vectors, and the transformation matrix as T(d_x, d_y, d_z) using the homogenous coordinates:

The matrix expression of the Equation (3.7) will be P' = P + T:

Using Equation (3.8), the location of the physical objects after any translation can be calculated easily. The format of the MDL function TRANSLATE is defined in the Table 3.9-1.

Table 3.9-1 Definition of The Function TRANSLATE

In a MCell simulation, the physical object may be rotated through an angle ø around any arbitrary axis (a, b, c). The axis (a, b, c) is the line extended from the world origin (0, 0, 0) to a given point (a, b, c). Positive angles of rotation are measured counter clockwise from the X axis toward the Y axis (right-hand rule). Derivation of the rotation matrix R(ø) of the rotation around the X, or Y, or Z axis can be found in computer graphics books and mathematical handbooks. Rotations around any arbitrary axis may be computed from the composition of the rotations around the X, Y, Z axis. The following is a short introduction about rotating around any arbitrary axis.

Let's denote point P(x, y, z) as a point on the selected physical object. After rotation through angle ø, it will become point P'(x', y', z'). If the rotation is around the X axis, the rotation matrix is R_x(ø):

similarly the rotation around the Y axis is expressed by the matrix R_y(ø), which is:

The rotation around the Z axis is defined by the matrix R_z(ø):

A rotation around any axis is computed from the combination of the matrices R_x, R_y, and R_z through five steps in the simulation. For a rotation of angle ø through an axis [a, b, c], we define a point A which is located at A(a, b, c), then the line segment OA is the rotation axis. The five steps of the computation about this arbitrary rotation are:

After the sequence of rotations, the point P(x, y, z) in the selected object is transformed to the point P'(x', y', z'), the coordinates of the point P'(x', y', z') can be computed from the equation as:

Any rotation around an axis [a, b, c] through an angle ø is defined by the MDL function ROTATE, the format of this function is shown in Table 3.9-2.

Table 3.9-2 Definition of the Function ROTATE

Beside translation and rotation, you can also scale the physical object in three dimensions. If the point P(x, y, z) is a point on the selected object, then after the object is scaled by a factor s_x along the X direction, a factor s_y along the Y direction, and a factor s_z along the Z direction, point P will become the point P'(x', y', z'). The scaling operation helps you to enlarge or shrink the object in the desired direction. Scaling is expressed by the equation P' = S(s_x, s_y, s_z)P. Where S is the scaling matrix with the format:

To have a matrix form of the scaling equation, we express the point P(x, y, z) and P'(x', y', z') in a vector format using the homogenous coordinates:

Then the matrix form of the scaling equation is:

The scaling shown in equation 3.18 is an operation about the origin of the Cartesian coordinate system. The physical object may be scaled by applying equation 3.16 to each point of the object, the proportion of the scaled object depends on the scaling factors in each direction, with a uniform scaling (s_x = s_y = s_z), the proportions of the scaled object will not be changed. Successive scaling are multiplicative, if you scale an object twice by the scaling matrices S1(s1_x, s1_y, s1_z), S2(s2_x, s2_y, s2_z), when the scaling operations are finished, the point P on the selected object will be turned to the point P''(x'', y'', z''), the point P'' can be computed from the equation:

The scaling operation is defined by the MDL command SCALE, the format of the command is given in Table 3.9-3.

Table 3.9-3 Definition of The Function SCALE

In order to get simulation results, you must specify the output format of the simulation in your input file. MCell can output the simulation results for a variety of analytical purposes, visualization data for graphical animation, numerical data for chemical reaction analysis, or the C language style output for the data analysis in a customized way.

The user may trace the statistics of the simulation event in a space or time dependent way. The amount of the simulation results are enormous, but you can specify the directory and the prefix of the output files in the MDL output specification.

MCell is able to export simulation results to fit the native data formats of a variety of graphic tools: Data Explore, IRIT, Rayshade, Renderman and Povray. To define the visualization type of results from the MCell simulation, you should choose your preferred visualization mode first, then specify when you want to save your results, and which object you want to observe. The MDL function VIZ_DATA_OUTPUT is used to define the output of the visualization results from the simulation. The format of the function VIZ_DATA_OUTPUT is defined in Table 3.10-1.

Table 3.10-1 The MCell Specification of The Visualization Output

You may use either the iteration number or the real simulation time to specify the simulation output of a MCell simulation. The two different output list formats used in MCell are: iteration_specifier, which uses iteration numbers to select the output range and time_specifier, which uses real time to define the output range. By using time_specifier you do not need to recalculate the iteration numbers for output when the time step of simulation changes. Inside a VIZ_DATA_OUTPUT function, you can only use one of the output_specifier, either the iteration_specifier or the time_specifier, you are not allowed use them both at the same time.

Iteration specifier (iteration_specifier) includes two MDL specifications that are used to select the iteration steps for output. The two iteration_specifier are ITERATION_LIST and ITERATION_FRAME_DATA. When you choose ITERATIONS_LIST as the iteration_specifier, all the information (ligand, effector, surface elements, positions, states etc.) of the selected objects will be saved at all the selected steps. On the contrary, the function ITERATION_FRAME_DATA is more flexible, it helps you to limit your output by specifying the interested properties of the objects at the specified iteration steps, so you will only get the results that you defined within the frame_data_specifier.

iteration_specifier lists the iteration numbers specified by the user for output. If we let the interger number N represent the iteration number, then the iteration_specifier has the format:

1. [N₁, N₂, N₃, ..., N_n]

2. [[N₁ TO N_n STEP dt]]

3. [list_of_iterations,[N₁ TO N_n STEP dt]]

4. [[N₁ TO N_n STEP dt], list_of_iterations]

The value of the MDL command ITERATIONS_LIST is a numerical array of the iteration times or an iteration_range_specifier, or the mixture of the iteration number lists and iteration_range_specifier. An iteration_range_specifier has the form as:

The TIME_LIST has exactly the same format as ITERATION_LIST, but instead of the iteration numbers it uses time to define the output. You may use a list of times, or a time range, or the combination of list and range of time to choose the simulation output.

time_specifier has the following possible values:

1. [t₁, t₂, t₃, ..., t_n]

2. [[t₁ TO t_n STEP dt]]

3. [list_of_times,[t₁ TO t_n STEP dt]]

4. [[t₁ TO t_n STEP dt], list_of_times]

Where t₁, t₂, t₃, t_n and dt are real numbers, they represent the time in the real world with a unit of second.

ITERATION_FRAME_DATA {

frame_data_specifier = iteration_specifier

... }

The MDL function ITERATION_FRAME_DATA specifies the type of the visualization data to be output and the simulation iteration numbers after which visualization data will be output. By associating the specified visualization data with indicated iteration numbers, the user can specify as many frame data associations as needed. The frame_data_specifierincludes:

The function OBJECT_FILE_PREFIXES is the file designator that is associated with the selected physical objects. This function selects the named surface objects for visualization output. The format of the function OBJECT_FILE_PREFIXES is:

OBJECT_FILE_PREFIXES {

surface_object_name = "output_domain_name"

...

}

In the definition of the function OBJECT_FILE_PREFIXES, output_domain_name is the user defined domain name that will be used to save the output files. You may define as many physical objects as you need in this function by giving them a different output_domain_name. The surface object name should include both the instantiated object name and the object name with a format like "instantiated_object_name.surface_object_name". Large amounts of data will be output to different files under the indicated directory or the current directory. If we use variable n to represent the selected iteration numbers in visualization output, then we will have the following output files:

For each iteration MCell generates a set of files like the one listed above.

Table 3.10-2 The MCell Chemical Reaction Output

The MDL command COUNT can be added as needed, you can combine the counted items to freely define your chemical reaction output. Now, let us see what we can count, where can we count and how can we count in the COUNT statements.

The reaction_specifier has three different expressions. It could be either iteration_specifier, time_specifier or STEP expression. We introduced iteration_specifier and time_specifier in the last section with output_specifier. Both the STEP and TIME_LIST definitions use time intervals (seconds) to count the points at which data is saved during the simulation iterations.

t, a real number, represents the time interval between each data output. This function lets the user save the reaction output with a fixed frequency, counting from the beginning of the simulation.

specifies the name of the reaction intermediate (ligand or effector state) or reaction state transition ( binding, unbinding, etc.), or the collision between ligand and surface to be output.

specifies the spatial location of the item indicated by count_what , it might be WORLD (the whole simulation world), existing_object or existing_region.

can be a list of specialized options for additional control over the manner of counting.

The allowed values for count_where and count_how depend upon the item specified by count_what. Note that mathematical expressions involving multiple COUNT statements may be specified within the braces {} surrounding the COUNT statement block. Please see the presently implemented COUNT statement blocks (subject to change and expansion):

The MCell simulator supports several generic file output commands which conform to the normal C language syntax. Please notice that only the file output statements are supported by MCell. We only list some commonly used output file functions in the following list. For more detailed description about the C language output functions please see the C language documentation.