14.1.4.5. Fick’s Law Data Model Scripting Example

MCell and CellBlender can be used to demonstrate various aspects of Fick’s 1st and 2nd laws. The Advanced Models section contains a tutorial on constructing such a model in CellBlender by hand. The following Data Model Script constructs an entire Fick’s demonstration model in Python.

As shown in the animation below, the script constructs a long semi-transparent blue box which will contain the molecules in the simulation. It then constructs a series of tan counting boxes along the length to count the molecules. It also constructs a series of light blue counting planes to measure the flux. Molecules are released at the center of the box at t=0 and diffuse during the simulation.

../_images/ficks_scripted.gif

14.1.4.5.1. Fick’s Law Construction Script

The following script builds the CellBlender model from scratch.

# Import cellblender to get the data model and project directory functions

import cellblender as cb

# Get a copy of the data model and change directories to the project directory for any file I/O

dm = cb.get_data_model(geometry=True)
old_location = cb.cd_to_project()

# Print some information about the data model

print ( "###############################################################" )

print ( "Data Model Top Level Keys = " + str(dm.keys()) )
print ( "MCell Keys = " + str(dm['mcell'].keys()) )

print ( "###############################################################" )


# Define parameters and set default values (these will also be available as local variables in this script)
#  The format for each parameter is:   Name, Value, Units, Description
#    All values are given as strings!!!

pars = [
  # These control the geometry of the model
  ['n',    "40",    '',     'Number of segments to sample along the x dimension'],
  ['lx',   "2.0",   'um',   'Length x = Total length of sample box in x dimension'],
  ['ly',   "0.2",   'um',   'Length y = Total length of sample box in y dimension'],
  ['lz',   "0.2",   'um',   'Length z = Total length of sample box in z dimension'],
  ['ext',  "0.02",  'um',   'Extended length for counting boxes and counting planes'],
  ['tol',  "0.995", 'um',   'Scale factor of sampling boxes to avoid coincident faces (0.995 works well)'],
  ['rtol', "0.001", 'um',   'Release plane tolerance - smaller is closer to ideal (0.001 works well)'],
  # These control the instrumentation of the model
  ['plot_segment_counts',   "1",  '',   'Plot count of vm molecules in each segment when non-zero'],
  ['plot_front_crossings',  "0",  '',   'Plot front crossings of vm molecules when non-zero'],
  # These describe the behavior of the model
  ['dc',    "5e-6",   'cm^2/sec',    'Diffusion Constant of Molecules'],
  ['nrel',  "1000",   'Count',       'Number of molecules to release'],
  # These control the simulation itself
  ['iters', "600",   '',      'Number of iterations to run'],
  ['seeds', "10",    '',      'Number of seeds to run'],
  ['dt',    "1e-6",  'sec',   'Time step for each iteration of the simulation'],
 ]

# Update local parameter list values from the existing data model with user-modified settings BEFORE regenerating it

dm_par_list = dm['mcell']['parameter_system']['model_parameters']
for dm_par in dm_par_list:
    print ( "Data Model Parameter " + dm_par['par_name'] + " = " + dm_par['par_expression'] )
    for p in pars:
        if dm_par['par_name'] == p[0]:
            # Update the local expression based on the parameter found in the incoming data model
            p[1] = dm_par['par_expression']
            p[2] = dm_par['par_units']
            p[3] = dm_par['par_description']


# Create the local variables from the updated values to use in this script
for p in pars:
    locals()[p[0]] = eval(p[1])


# Create the new mcell data model inside the existing data model (this deletes the previous mcell data model)

dm['mcell'] = { 'data_model_version' : "DM_2014_10_24_1638" }


# Restore the parameters that were either initialized from scratch or preserved from the previous data model

dm['mcell']['parameter_system'] = { 'model_parameters':[] }   # Parameters are currently unversioned
for p in pars:
    dm['mcell']['parameter_system']['model_parameters'].append ( { 'par_name':p[0], 'par_expression':p[1], 'par_units':p[2], 'par_description':p[3] } )


# Define a function to make either a plane or a box from its center and lengths in each dimension (one zero dimension gives a plane)

def make_obj ( center_x, center_y, center_z, len_x, len_y, len_z ):
    obj = {}
    obj['vertex_list'] = []
    obj['element_connections'] = []

    if len_x == 0:
      # Make a plane perpendicular to the x axis
      obj['vertex_list'].append ( [ center_x, center_y-(len_y/2.0), center_z-(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x, center_y-(len_y/2.0), center_z+(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x, center_y+(len_y/2.0), center_z+(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x, center_y+(len_y/2.0), center_z-(len_z/2.0) ] )
      obj['element_connections'].append ( [ 0, 2, 1 ] )
      obj['element_connections'].append ( [ 0, 3, 2 ] )
    elif len_y == 0:
      # Make a plane perpendicular to the y axis
      obj['vertex_list'].append ( [ center_x-(len_x/2.0), center_y, center_z-(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x-(len_x/2.0), center_y, center_z+(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x+(len_x/2.0), center_y, center_z+(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x+(len_x/2.0), center_y, center_z-(len_z/2.0) ] )
      obj['element_connections'].append ( [ 0, 2, 1 ] )
      obj['element_connections'].append ( [ 0, 3, 2 ] )
    elif len_z == 0:
      # Make a plane perpendicular to the z axis
      obj['vertex_list'].append ( [ center_x-(len_x/2.0), center_y-(len_y/2.0), center_z ] )
      obj['vertex_list'].append ( [ center_x-(len_x/2.0), center_y+(len_y/2.0), center_z ] )
      obj['vertex_list'].append ( [ center_x+(len_x/2.0), center_y+(len_y/2.0), center_z ] )
      obj['vertex_list'].append ( [ center_x+(len_x/2.0), center_y-(len_y/2.0), center_z ] )
      obj['element_connections'].append ( [ 0, 2, 1 ] )
      obj['element_connections'].append ( [ 0, 3, 2 ] )
    else:
      # Make a box
      obj['vertex_list'].append ( [ center_x+(len_x/2.0), center_y+(len_y/2.0), center_z-(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x+(len_x/2.0), center_y-(len_y/2.0), center_z-(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x-(len_x/2.0), center_y-(len_y/2.0), center_z-(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x-(len_x/2.0), center_y+(len_y/2.0), center_z-(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x+(len_x/2.0), center_y+(len_y/2.0), center_z+(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x+(len_x/2.0), center_y-(len_y/2.0), center_z+(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x-(len_x/2.0), center_y-(len_y/2.0), center_z+(len_z/2.0) ] )
      obj['vertex_list'].append ( [ center_x-(len_x/2.0), center_y+(len_y/2.0), center_z+(len_z/2.0) ] )
      obj['element_connections'].append ( [ 1, 2, 3 ] )
      obj['element_connections'].append ( [ 7, 6, 5 ] )
      obj['element_connections'].append ( [ 4, 5, 1 ] ) # Right end
      obj['element_connections'].append ( [ 5, 6, 2 ] )
      obj['element_connections'].append ( [ 2, 6, 7 ] ) # Left end
      obj['element_connections'].append ( [ 0, 3, 7 ] )
      obj['element_connections'].append ( [ 0, 1, 3 ] )
      obj['element_connections'].append ( [ 4, 7, 5 ] )
      obj['element_connections'].append ( [ 0, 4, 1 ] ) # Right end
      obj['element_connections'].append ( [ 1, 5, 2 ] )
      obj['element_connections'].append ( [ 3, 2, 7 ] ) # Left end
      obj['element_connections'].append ( [ 4, 0, 7 ] )

    return obj


# Add materials for the objects

dm['mcell']['materials'] = { 'material_dict' : {} }   # Materials are currently unversioned
dm['mcell']['materials']['material_dict']['box_color']   = { 'diffuse_color' : {'a':0.3, 'r':0.2, 'g':0.4, 'b':1.0} }
dm['mcell']['materials']['material_dict']['rel_color']   = { 'diffuse_color' : {'a':0.2, 'r':0.9, 'g':0.7, 'b':0.5} }
dm['mcell']['materials']['material_dict']['vol_color']   = { 'diffuse_color' : {'a':0.1, 'r':0.9, 'g':0.7, 'b':0.5} }
dm['mcell']['materials']['material_dict']['plane_color'] = { 'diffuse_color' : {'a':0.7, 'r':0.5, 'g':0.7, 'b':1.0} }



# Create container objects for geometrical objects and model objects

dm['mcell']['geometrical_objects'] = {}   # Geometrical objects are currently unversioned
dm['mcell']['model_objects'] = { 'data_model_version':"DM_2014_10_24_1638" }

# Each container also includes a list

dm['mcell']['geometrical_objects']['object_list'] = []
dm['mcell']['model_objects']['model_object_list'] = []


# Add objects to the lists

# Make the main box for diffusing the molecules

box = make_obj ( 0, 0, 0, 10*lx, ly, lz ) # Make the box much longer to reduce boundary effects from absorptive ends
box['name'] = 'box'
box['material_names'] = [ 'box_color' ]

# Make the thin box for releasing the molecules

rel = make_obj ( 0, 0, 0, rtol, ly-rtol, lz-rtol )
rel['name'] = 'rel'
rel['material_names'] = [ 'rel_color' ]

# Make the surface regions for the two absorptive ends to keep them from accumulating
box['define_surface_regions'] = []
box['define_surface_regions'].append ( { 'name':"left_end", 'include_elements':[ 4, 10 ] } )
box['define_surface_regions'].append ( { 'name':"right_end", 'include_elements':[ 2, 8 ] } )

# Add the box to the geometrical objects and the model objects
dm['mcell']['geometrical_objects']['object_list'].append ( box )
dm['mcell']['model_objects']['model_object_list'].append ( { 'name':box['name'] } )

dm['mcell']['geometrical_objects']['object_list'].append ( rel )
dm['mcell']['model_objects']['model_object_list'].append ( { 'name':rel['name'] } )

# Make the counting boxes and planes as requested by the parameter flags
for i in range(n):

    x = (i - ((n-1)/2.0)) / (n/lx)

    if plot_segment_counts != 0:
      box = make_obj ( x, 0, 0, tol*(lx/n), ly+ext, lz+ext )
      box['name'] = 'vol_%03d' % i
      box['material_names'] = [ 'vol_color' ]
      dm['mcell']['geometrical_objects']['object_list'].append ( box )
      dm['mcell']['model_objects']['model_object_list'].append ( { 'name':box['name'] } )

    if (plot_front_crossings != 0) and (i > 0):
      plane = make_obj ( x-(lx/(2*n)), 0, 0, 0.0, ly+ext+ext, lz+ext+ext )
      plane['name'] = 'plane_%03d' % i
      plane['material_names'] = [ 'plane_color' ]
      dm['mcell']['geometrical_objects']['object_list'].append ( plane )
      dm['mcell']['model_objects']['model_object_list'].append ( { 'name':plane['name'] } )


# Create a molecule list and create a "vm" molecule along with its display properties in that list

dm['mcell']['define_molecules'] = { 'data_model_version' : "DM_2014_10_24_1638" }
mol = { 'mol_name':"vm", 'mol_type':"3D", 'diffusion_constant':"dc", 'data_model_version':"DM_2016_01_13_1930" }
mol['display'] = {'color':[0.0,1.0,0.0], 'emit':1.0, 'glyph':"Cube", 'scale':0.5 }
dm['mcell']['define_molecules']['molecule_list'] = [ mol ]


# Create a release site

rel_site = {
              'name' : "center_rel",
              'molecule' : "vm",
              'quantity' : "nrel",
              'quantity_type' : "NUMBER_TO_RELEASE",
              'release_probability' : "1",
              'shape' : "OBJECT",
              'object_expr' : "rel",
              'orient' : ";",
              'pattern' : "",
              'location_x' : "0",
              'location_y' : "0",
              'location_z' : "0",
              'site_diameter' : "0",
              'stddev' : "0",
              'data_model_version' : "DM_2015_11_11_1717"
           }

dm['mcell']['release_sites'] = { 'release_site_list':[ rel_site ], 'data_model_version':"DM_2014_10_24_1638" }

# Define surface classes

dm['mcell']['define_surface_classes'] = { 'surface_class_list':[], 'data_model_version':"DM_2014_10_24_1638" }

# Use a table to construct the various classes with associated properties

surf_classes = [
  [ 'transp', 'vm_transp', ';', "TRANSPARENT",         "0" ],
  [ 'absorb', 'vm_absorb', ';', "ABSORPTIVE",          "0" ] ]

# Loop through the table and add each class to the data model

for c in surf_classes:
  sc_prop = { 'data_model_version':"DM_2015_11_08_1756",
              'name':c[1],
              'affected_mols':"SINGLE",
              'molecule':"vm",
              'surf_class_orient':c[2],
              'surf_class_type':c[3],
              'clamp_value':c[4]
            }

  sc_entry = { 'data_model_version':"DM_2014_10_24_1638",
               'name':c[0],
               'surface_class_prop_list':[ sc_prop ]
             }

  dm['mcell']['define_surface_classes']['surface_class_list'].append ( sc_entry )


# Assign the surface classes with the "modify_surface_regions" key

dm['mcell']['modify_surface_regions'] = { 'modify_surface_regions_list':[], 'data_model_version': "DM_2014_10_24_1638" }


# Modify the left end to be absorptive

dm['mcell']['modify_surface_regions']['modify_surface_regions_list'].append (
    {
      'name':"absorb left",
      'object_name':"box",
      'region_name':"left_end",
      'surf_class_name':"absorb",
      'region_selection':"SEL",
      'data_model_version':"DM_2015_11_06_1732"
    } )

# Modify the right end to be absorptive

dm['mcell']['modify_surface_regions']['modify_surface_regions_list'].append (
    {
      'name':"absorb right",
      'object_name':"box",
      'region_name':"right_end",
      'surf_class_name':"absorb",
      'region_selection':"SEL",
      'data_model_version':"DM_2015_11_06_1732" } )

# Modify the release box, all counting boxes, and counting planes (if any) to be transparent

dm['mcell']['modify_surface_regions']['modify_surface_regions_list'].append (
  {
    'name':"transp rel",
    'object_name':"rel",
    'region_name':"",
    'surf_class_name':"transp",
    'region_selection':"ALL",
    'data_model_version':"DM_2015_11_06_1732"
  } )

for i in range(n):

    if plot_segment_counts != 0:
      name = 'vol_%03d' % i
      dm['mcell']['modify_surface_regions']['modify_surface_regions_list'].append (
        {
          'name':"transp "+name,
          'object_name':name,
          'region_name':"",
          'surf_class_name':"transp",
          'region_selection':"ALL",
          'data_model_version':"DM_2015_11_06_1732"
        } )

    if (plot_front_crossings != 0) and (i > 0):
      name = 'plane_%03d' % i
      dm['mcell']['modify_surface_regions']['modify_surface_regions_list'].append (
        {
          'name':"transp "+name,
          'object_name':name,
          'region_name':"",
          'surf_class_name':"transp",
          'region_selection':"ALL",
          'data_model_version':"DM_2015_11_06_1732"
        } )


# Define the counting output

dm['mcell']['reaction_data_output'] = {
    'data_model_version':"DM_2014_10_24_1638",
    'reaction_output_list':[],
    'rxn_step':"10*dt",
    'combine_seeds':False,
    'mol_colors':True,
    'plot_layout':" plot ",
    'plot_legend':"x",
    'mol_colors':False
    }

dm['mcell']['reaction_data_output']['reaction_output_list'].append (
    {
      'data_model_version':"DM_2015_10_07_1500",
       'name':"vm in box",
       'rxn_or_mol':"Molecule",
       'mdl_string':"",
       'mdl_file_prefix':"",
       'count_location':"Object",
       'object_name':"box",
       'region_name':"",
       'reaction_name':"",
       'molecule_name':"vm"
    } )

# Create the counting structures for the counting object as requested

for i in range(n):
    if plot_segment_counts != 0:
        name = 'vol_%03d' % i
        if plot_segment_counts != 0:
          dm['mcell']['reaction_data_output']['reaction_output_list'].append (
            {
              'data_model_version':"DM_2015_10_07_1500",
              'name':"vm in "+name,
              'rxn_or_mol':"Molecule",
              'mdl_string':"",
              'mdl_file_prefix':"",
              'count_location':"Object",
              'object_name':name,
              'region_name':"",
              'reaction_name':"",
              'molecule_name':"vm"
            } )
    if (plot_front_crossings != 0) and (i > 0):
      name = 'plane_%03d' % i
      mdl_string = "COUNT[vm,Scene."+name+",FRONT_CROSSINGS]"
      dm['mcell']['reaction_data_output']['reaction_output_list'].append (
        {
          'data_model_version':"DM_2015_10_07_1500",
          'name':"MDL: "+mdl_string,
          'rxn_or_mol':"MDLString",
          'mdl_file_prefix':name+"_front_cross",
          'mdl_string':mdl_string,
          'count_location':"World",
          'object_name':"",
          'region_name':"",
          'reaction_name':"",
          'molecule_name':""
        } )


# Set up the simulation running parameters

dm['mcell']['initialization'] = { 'data_model_version':"DM_2014_10_24_1638" }
dm['mcell']['initialization']['iterations'] = "iters"
dm['mcell']['initialization']['time_step'] = "dt"

dm['mcell']['simulation_control'] = { 'data_model_version': 'DM_2016_04_15_1430' }
dm['mcell']['simulation_control']['start_seed'] = '1'
dm['mcell']['simulation_control']['end_seed'] = 'seeds'


# Return to the previous directory and replace the existing data model with this modified version

cb.cd_to_location ( old_location )
cb.replace_data_model ( dm, geometry=True )

14.1.4.5.2. Fick’s Law Plotting Script

The following plot series (animation) compares the average of 50 MCell runs (Start Seed=1, End Seed=50, 10,000 molecules released) at various points in time to plots of the expected theoretical results at those same points in time.

../_images/plot_ideal_mcell_both_crop.gif

This script generates the data needed for that plot. It reads and averages the data from the MCell runs and generates a series of plot files from that output. It also generates a series of plot files from the expected analytic solution with the same parameters as the MCell runs. Because this script uses the same parameters used to generate the simulation, it will produce the correct plot even when the simulation parameters are varied.

import math
import os
from numpy import fromfile

import cellblender as cb

dm = cb.get_data_model()
old_location = cb.cd_to_project()

pars = [
  # These control the geometry of the model
  ['n',    "40",    '',     'Number of segments to sample along the x dimension'],
  ['lx',   "2.0",   'um',   'Length x = Total length of sample box in x dimension'],
  ['ly',   "0.2",   'um',   'Length y = Total length of sample box in y dimension'],
  ['lz',   "0.2",   'um',   'Length z = Total length of sample box in z dimension'],
  ['ext',  "0.02",  'um',   'Extended length for counting boxes and counting planes'],
  ['tol',  "0.995", 'um',   'Scale factor of sampling boxes to avoid coincident faces (0.995 works well)'],
  ['rtol', "0.001", 'um',   'Release plane tolerance - smaller is closer to ideal (0.001 works well)'],
  # These control the instrumentation of the model
  ['plot_segment_counts',   "1",  '',   'Plot count of vm molecules in each segment when non-zero'],
  ['plot_front_crossings',  "0",  '',   'Plot front crossings of vm molecules when non-zero'],
  # These describe the behavior of the model
  ['dc',    "5e-6",   'cm^2/sec',    'Diffusion Constant of Molecules'],
  ['nrel',  "1000",   'Count',       'Number of molecules to release'],
  # These control the simulation itself
  ['iters', "500",   '',      'Number of iterations to run'],
  ['seeds', "10",    '',      'Number of seeds to run'],
  ['dt',    "1e-6",  'sec',   'Time step for each iteration of the simulation'],
 ]

# Update local parameter list values from the existing data model with user-modified settings BEFORE generating plot data

dm_par_list = dm['mcell']['parameter_system']['model_parameters']
for dm_par in dm_par_list:
    print ( "Data Model Parameter " + dm_par['par_name'] + " = " + dm_par['par_expression'] )
    for p in pars:
        if dm_par['par_name'] == p[0]:
            # Update the local expression based on the parameter found in the incoming data model
            p[1] = dm_par['par_expression']
            p[2] = dm_par['par_units']
            p[3] = dm_par['par_description']

# Create the local variables from the updated values to use in this script
for p in pars:
    locals()[p[0]] = eval(p[1])


plot_iters = [ 10, 25, 50, 100, 200, 300, 400, 600 ]
start_seed = 1
end_seed = seeds
start_vol = 0
end_vol = n - 1
num_vols = 1 + end_vol - start_vol

# Note: Neither 'iterations' nor 'time_step' can use CellBlender parameters!!
iters = eval(dm['mcell']['initialization']['iterations'])
dt = eval(dm['mcell']['initialization']['time_step'])

react_files_dir = "mcell" + os.sep + "react_data"
react_files_seeds = os.listdir(react_files_dir)

print ( "files: " + str(react_files_seeds) )

parent_dir = os.getcwd().split(os.sep)[-1]

# Remove all plots from the data model so they're not cumulative when this script is re-run

dm['mcell']['reaction_data_output']['reaction_output_list'] = []  # Comment this line to NOT remove previous plots


def make_file_plot ( file_name ):
  blender_relative_name = "//" + parent_dir + os.sep + fname
  reaction_out = { 'data_model_version':"DM_2016_03_15_1800" }
  reaction_out['rxn_or_mol'] = "File"
  reaction_out['molecule_name'] = ""
  reaction_out['reaction_name'] = ""
  reaction_out['object_name'] = ""
  reaction_out['region_name'] = ""
  reaction_out['mdl_file_prefix'] = ""
  reaction_out['mdl_string'] = ""
  reaction_out['count_location'] = "World"
  reaction_out['plotting_enabled'] = True
  reaction_out['data_file_name'] = blender_relative_name
  reaction_out['name'] = "FILE:" + blender_relative_name
  return reaction_out


area = ly * lz

# Generate the MCell plots

for plot_iter in plot_iters:
  print ( "Generating concentration curve for " + str(plot_iter) )

  fname = "concentration_%d.txt" % plot_iter
  f = open ( fname, "w" )
  dm['mcell']['reaction_data_output']['reaction_output_list'].append ( make_file_plot(fname) )

  points = []
  for vol in range(start_vol,1+end_vol):

      count = 0.0
      sx0 = (1+end_vol-start_vol)/2
      sx = (vol - sx0)
      vx = sx * lx / n
      for seed in range(start_seed,1+end_seed):
          file_name = react_files_dir + os.sep + ("seed_%05d" % seed) + os.sep + ("vm.vol_%03d.dat" % vol)
          data = fromfile ( file_name, sep=' ' )
          x = data[0::2]
          y = data[1::2]
          count = count + y[plot_iter]
      averaged_count = count/(1+end_seed-start_seed)
      conc = averaged_count / (tol*lx/n)                    # This should be molecules per micron (length)

      f.write ( str(vx+(lx/(2.0*n))) + "  " + str(conc) + "\n" )
  f.close()


# Generate the analytic plots

for plot_iter in plot_iters:
  print ( "Generating analytic curve for " + str(plot_iter) )

  fname = "concentration_ideal_%d.txt" % plot_iter
  f = open ( fname, "w" )
  dm['mcell']['reaction_data_output']['reaction_output_list'].append ( make_file_plot(fname) )

  t = plot_iter * dt

  for vol in range(start_vol,1+end_vol):
      sx0 = (1+end_vol-start_vol)/2
      sx = (vol - sx0)
      x = sx * lx / n
      N = nrel * math.exp(-(x*x/(4*dc*1e8*t))) / (2 * math.sqrt(math.pi*dc*1e8*t))
      f.write ( str(x) + "  " + str(N) + "\n" )
  f.close()

# Configure the plotter for one page / one plot (uses a single space):

dm['mcell']['reaction_data_output']['plot_layout'] = ' '

cb.cd_to_location ( old_location )
cb.replace_data_model ( dm )