Modeling Overview

This document briefly explains the modeling process, from microscope slide to experimental results

Geometric Model

Original instructions

It is probably ok to delete these. New instructions follow.

=====================================================
To use the mesh, you need a series of Matlab files.

nerve_model.m
 this is a script, so just call it, and it is as if you are typing the
everything at the command line. When it is done, you have read in the
bitmap, flattened it, and replaced RGB values with sequential integers.

make_myelin_table.m
 this is a function, so look at the top line of the code to see what
variables it needs. I named them such that the variables in the
function are the same as the ones made at the end of the nerve_model
script. Therefore, you could cut and paste the function header if you
wanted to, and end up with the right arguments. The output is
myelin_table.

make_axon_table.m
 same as above, in terms of structure and use. Output is axon_table.

extrude.m
 similar to above, but the last argument in the function is not a
variable; instead, you need to put in how thick you want the extrusion
to be. The output is extruded_model.m

you can linearize the model with the linearize function in the
linearize.m file, if you want to.

Matrix conventions:
Matlab matrices are ordered M(row,column). Our digitized bitmaps are oriented x horizontal, y vertical, consistent with Gimp and Photoshop. So the correct order when transferring data is M(y,x).

>> nerve_model
>> myelin_table=make_myelin_table(indices,nerve_map);
>> axon_table=make_axon_table(myelin_table)
>> extruded_model=extrude(nerve_map,axon_table,myelin_table,indices,10)

NOTE: you can substitute any integer for 10; it is how thick the model
will be.

=====================================================

New instructions, June 19, 2003

% change the directory to where the working scripts are located

cd matlab_files/working_scripts/

% run the “nerve_model” script, which is a .m file, but is not a

% function like the other .m files are. It is a script file, so each command

% is treated as if it were entered at the command line.

% usually takes ~5 minutes

nerve_model

% use the indices and nerve_map matrices that are generated above to make a

% table that catalogs all of the myelin in the picture. myelin_table will be

% a record x and y locations for every myelin pixel that is in the original

% bitmap, and axoplasm table will be a record x and y locations for every

% axoplasm pixel that is in the original

% usually takes ~5 minutes.

myelin_table=make_myelin_table(indices,nerve_map);

axoplasm_table=make_axoplasm_table(indices,nerve_map);

% aggregate the data in the “myelin_table” matrix and store in the

% “myelinated_axon_table”, so that all of the *contiguous* myelin is assumed to

% be part of the same axon, so that mean x and y values, diameters, and other

% descriptive characteristics can be generated for each axon.

% do the same for the “axoplasm_table” matrix, and store the result in the

% “all_axon_table”.

% usually takes ~15 minutes

myelinated_axon_table=make_axon_table(myelin_table);

all_axon_table=make_axon_table(axoplasm_table);

% make the transfer function tables, which are the lensed and reduced m and n

% values for each x and y position. By definition, any lensed coordinate is an

% <m,n> pair, which can be either reduced or unreduced

% The warning regarding corner preservation simply refers to what happens to

% pixels that aren’t within the area of the lens – if the corners are preserved,

% they are just transferred directly to the lensed bitmap. If they are not, the

% corners are filled with 0 values. In general, preservation is what is

% desired.

transfer_table=transfer_func(nerve_map);

% make the invers_transfer_table, which will look up the unreduced <x,y> pair for each

% reduced <m,n> pair.

inverse_transfer_table=invert_transfer_table(transfer_table);

% extrude the model in the z direction, store the z values in the z_coordinates matrix,

% and update the indices array

% New output added by crb, September 2003: node_list is a list of xyz centers of nodes of Ranvier

[extruded_model,z_coordinates,indices,node_list]=extrude(nerve_map, transfer_table, inverse_transfer_table, myelinated_axon_table, myelin_table, indices);

% Run register.m to create a lookup table from all_axon_table to myelinated_axon_table

% Added 8/1/03, crb

registration_key=register(all_axon_table,myelinated_axon_table);

% now correct the seed values in all_axon_table so that myelinated fibers values are correct

% Added 8/1/03, crb

all_axon_table(registration_key(:,1),9)=myelinated_axon_table(registration_key(:,2),9);

% ground the model by putting a new compartment in the corners, upon which boundary

% conditions can be imposed.

[grounded_extruded_model,indices]=ground_corners(extruded_model,indices);

% linearize the model, which will be un-linearized during importation into

% SCIRun, if the model geometry is correctly described during the importation

% process on that end

linear_extruded_model=linearize(grounded_extruded_model);

% save the results to a file; the all_axon_table and myelinated_axon_table will

% be necessary, in particular, for the post-hoc analysis.

save('../geometry_data/050203/linear_extruded_model_050203','linear_extruded_model');

save('../geometry_data/050203/linear_z_050203.mat', 'z_coordinates')

% make a “node transfer table”. The ones we made earlier were for the cells. In order to

% identify node locations for sci-run, the same routines must be executed for a matrix

% that is one pixel larger in the x and y direction than the original bitmap (if the

% original number of cells is 821x749, the number of nodes is 822x750). The actual data

% in the matrix is irrelevant, because the only “input” to the transfer_func is the

% x and y locations.

node_map=zeros(822,750);

node_transfer_table=transfer_func(node_map);

% linearize the x and y coordinates, and save them to a file for SCIRun.

node_coordinates=reshape(node_transfer_table,[24600,2]);

x_coordinates=node_coordinates(:,1);

y_coordinates=node_coordinates(:,2);

clear node_coordinates

save('../geometry_data/050203/linear_x_050203.mat', 'x_coordinates')

save('../geometry_data/050203/linear_y_050203.mat', 'y_coordinates')

save ../geometry_data/050203/nerve_model_050302.mat

Data Conversion

At this point there are three Matlab variables saved in a file:

-nerve_model_(date): contains conductivity indices for each pixel
-linear_x_(date).mat, linear_y_(date).mat, linear_z_(date).mat: contain x, y and z coordinates for each pixel.

Now run structhex_converter.m.  This Matlab script uses the three input files just mentioned to generate two new files in a format suitable for the scirun converter:
-structdata.nodes: contains xyz coordinates of each node in the mesh.  Note that the number of nodes is the number of pixels+1 in each dimension. Note also that the xyz positions of each node are rescaled to .29 microns/pixel.
-structdata.cond: Indicates conductivity index number or tissue type.  Actual conductivity values are added later.

These two input files are converted to a SCIRun mesh with:
RawToStructHexVol structdata.nodes structdata.cond structhex.fld

Note that structhex_converter.m also creates pointcloud files for the ground nodes.  Specifically, it creates:
-corners: Contains corner nodes for all slices.  Convert to scirun format with: TextToPointCloudField corners corners.fld
-top_bottom_point_cloud: contains the top and bottom slices.  This file probably won't be used any longer.  Instead, a network called  clip_ends.net is used to get the top and bottom from the outer mesh only.  This is done for two reasons: 1. it removes ground nodes from inside the axons and 2. it allows all stimulating nodes to be in the outer (non-axoplasm) mesh.

Biophysical Model

Once the mesh is in SCIRun format there are several steps to complete the biophysical model.

• Add conductivities: Use add_conds.net to add conductivity values.  As of 6/2/03 the conductivities in S/m are:
• 0: corner area. There aren’t any of these
• 1: unidentified tissue near the perimeter.  There are very few of these.
• 2: epineurium .1
• 3: endothelium .6
• 4: perineurium .01
• 5: intravascular .6
• 6: endoneurium .333
• 7: axoplasm .909
• 8: myelin 1.24e-6
• 9: nodes of Ranvier .333
• 10: corner ground nodes 1
• Prepare tetvols and surfaces with prep_fields.net, which does the following:
• Converts from StructHex to TetVol using the HexToTet module
• Saves the outer mesh using ClipByFunction v != 7 as tetvol_outer.fld
• Saves the inner mesh using ClipByFunction v == 7 as tetvol_axoplasm.fld
• Finds the axoplasm boundary by using FieldBoundary on the inner mesh.  Save as axoplasm_boundary.fld
• Prepare the cathode with clip_ends.net
• The cathode is made up of the top and bottom slices (minus the nodes within the axons) and the corner nodes.
• Clips the ends of the mesh with ClipByFunction z < -7000 | z > 7000 from tetvol_outer.fld. This must be performed twice: the first iteration produces a tetvol and the second iteration produces a pointcloud.
• Add the ends to the corners with GatherPoints and save as cathode.fld
• Prepare the stiffness matrices and mesh with prep_mesh_and_stiffmat.net
• SetupFEMatrix for inner and outer meshes, saved as axoplasm_stiffmat.mat and outer_stiffmat.mat respectively
• Adds cathodes to outer mesh, save as tetvol_outer_w_ground_electrodes.fld. Note that “Interpolate outside mesh” must be checked in InsertVoltageSource. Otherwise some of the corner nodes might get skipped.
• Build interpolant from axoplasm mesh to outer mesh with interp_outer_to_inner.net
• Save interpolant matrix as outer_to_inner_interp.mat
• Convert to text format with SparseRowMatrixToText outer_to_inner_interp.mat outer_to_inner_interp.txt -noHeader -oneBasedIndexing
• Calculate the membrane surface with clip_axoplasm_boundary.net, which does the following:
• Finds the cell surface between the outer and inner meshes.
• Saves the field data as axosurface_boundary.fld (this is different than the axoplasm volume surface, which includes the ends).
• Saves interpolant from outer mesh to axosurface and membrane surface area as Matlab file interp_outer_to_axosurface_and_face_area.m, where i1 is the interpolant and i2 is the membrane area vector.
• Convert the axosurface field with: TriSurfFieldToText axosurface_boundary_no_ends.fld axosurface_points axosurface_faces -noPtsCount -noElementsCount -oneBasedIndexing
• Prepare the membrane stiffness matrix quadrants with the Matlab script resistive_model.m which does the following:
• Creates 4 new matrices: top_left_stiffmat.txt, top_right_stiffmat.txt, bot_left_stiffmat.txt and bot_right_stiffmat.txt.
• Each of these must be converted to SCIRun format with TextToSparseRowMatrix (name)_stiffmat.txt (name)_stiffmat.mat -binOutput -oneBasedIndexing
• Creates bot_rhs.txt which is a blank right hand side representing sources in the inner mesh. It must be appended to the outer right hand side to solve the system.
• This file must be converted to SCIRun format with TextToSparseRowMatrix bot_rhs.txt bot_rhs.mat -binOutput
• Prepare the new stiffness matrix with membrane and solve the system using 3d_solution.net, which does the following:
• Add 4x7 array of anodes to outer mesh and injects current into one electrode. The anode array is 60x60x3000 microns. Note that ApplyFEMCurrentSource will only accept the current sources as a PointCloud; this is a problem because SampleLattice creates a LatVolField. While GatherPoints will convert the LatVol to a PointCloud, it strips out all of the data. So the temporary solution is to use GatherPoints with DirectInterpolate.
• Save anode positions to anode_xyz.mat for analysis purposes (distance from active electrode to fiber, node of Ranvier, etc.)
• Build new outer stiffness matrix with ground nodes and current sources.
• Build right hand side for outer mesh and expand it to accommodate inner mesh (i.e. increase its length)
• Assemble the final stiffness matrix by adding the new outer stiffness matrix to top_left, adding the inner to bot_right, and concatenating these results with top_right and bot_left.
• Run the solver to calculate the current density vector field. Use the following solver settings:
1. Emit partial solutions: off
2. Use previous solution as initial guess: off
3. Maximum iterations: >6000
4. Target error: 1e-4 (per Oleg, through Alexei)
• Save the vector positions in the current density vector field as vector_xyz.mat
• Save the magnitudes of the current density vector field as vector_ijk_anode<anode#>_<intensity value>uA.mat
• Note that the results must be converted to Matlab with the MatlabInterface – text converters will not work because the numbers are too small.
• Addendum, September 12, 2003
• In order to solve the system accurately, the petsc solvers must be used
• Petsc solvers only run on linux, single processor, which in the SCI lab means they must be run on inferno
• This means that the stiffness matrix and rhs must be prepared in one network, and the system must be solved in a second network because the solver takes all of the available memory (2 GB).
• The process so far has been to do the preparation with 3d_solver.net and solution with solve_matrix.net
• The plan at the moment is to solve for each electrode with 1uA current and scale the solution appropriately for other current values.
• Solver Notes
• The default conjugate gradient solver in SCIRun will not work, mainly because it performs an inadequate search of the eigenvector field. The consequence is that the results contain anomalous activation sites.
• Petsc solvers seem to offer much better results. There are about 15 solvers and 15 preconditioners available in SCIRun. I have tested about half of these combinations, and the PBCGS solver with the MILU preconditioner seems to provide the best results so far in the smallest number of iterations.
1. I have tried error levels of 1e-6, 1e-7 and 1e-8.
2. I think 1e-8 is the optimal error level.
• One additional subject is whether it is necessary to solve the system for each current level, since all results could be expressed in terms of the solution at 1uA. At the moment, it seems like only one solution per electrode is necessary.
1. I solved for electrode 13 with current injection levels of 1uA and 10uA. The results show that the z component of the current density vector field are 10 times larger for the latter. This is encouraging.

Notes on axon identification

• There is a problem with axon identification from x,y coordinates. It appears that the cell coordinates used throughout the analysis are slightly shifted from the hexvol center coordinates in SCIRun. But this is not a scirun problem – it results from the fact that the hexvol centers are determined from the node_map, and this results in slightly shifted coordinates. So all of the original analysis code to identify which neurons fired misidentified some neurons.
• The solution is to identify neurons based on pixel position rather than x,y coordinates, since the coordinates will not match exactly. This requires that the max current density in the z direction be mapped
• First calculate the current density vector field as explained above.

Notes on variable slice thickness in Z direction

• Variable slice thickness in the z direction changes the average size of the membrane faces, and thus changes the resistance (and capacitance) values for that patch of membrane. So, instead of using an average area size to calculate R and C, the area of each face must be determined. There are three approaches to solving this:
1. Clips ends off of the axoplasm boundary. This would prevent any membrane from being present at the ends of the model, which would not be a major problem but would create a mismatch between membrane nodes and intersection nodes.
2. Leave ends in the model, don't worry about nodes inside the axon. The error here would be small, but would also have a mismatch in nodes between membrane and intersection.
3. Remove center of axon on ends. This doesn't work for small axons that are 1 hexvol wide, but does work otherwise.
• The best of these options is 3. It has the smallest error, and more importantly preserves a 1 to 1 relationship between membrane nodes and intersection nodes.

Calculating distance from electrodes to fibers and nodes of Ranvier

electrode positions are stored in anode.txt (in scirun coordinates) so
load it first:

load anode.txt

node_list must be generated by extrude.m and is not stored in
nerve_model_060203.mat (or similar data file) so it may be necessary to
run extrude. Also, myelinated_axon_centers (in scirun coordinates) must
be loaded:
load myelinated_axon_centers

Now run:
[dist_to_fiber,dist_to_nor,node_xyz]=dist_to_electrode(anode,myelinated_axon_centers,node_list);

which returns all variables in scirun coordinates. Save the dist_to_fiber
and dist_to_nor tables:
save dist_to_fiber dist_to_fiber
save dist_to_nor dist_to_nor

Then import them into scirun and attach them to the proper fields. Three new fields have been added:

myelinated_axon_centers.fld

all_axon_centers.fld

all_axon_fiber_size.fld

Notes on Results

• It appears that it will only be necessary to solve for each electrode position once, for example at 1uA. The voltages at all other current levels are multiples of the 1uA solution.
• Further, it is only necessary to calculate the gradient of the field within the axoplasm. Since the conductivity of the axoplasm is homogenous and isotropic, calculating the current density vector field (gradient dot sigma) will just change the gradient by a scaling factor (-.909 in this case).

Notes on myelin properties

• Myelin has both resistance and capacitance, similar to the cell membrane. The major difference is that the cell membrane is extremely thin, while myelin has nonzero thickness (at least the way it is constructed in our model).
• We could insert myelin capacitance in a manner similar to membrane capacitance, except that the myelin has nonzero thickness and variable resistivity. So the capacitance is proportional not only to the surface area but also the thickness.
• McIntyre addressed the myelin problem by creating double cable models (using the cable equation).

Notes on nodes of Ranvier indices

· Nodes of Ranvier are spaced at 100*fiber diameter and then rounded to the nearest slice. Unfortunately, all_axon_table and myelinated_axon_table contain different seed values for nodes of Ranvier. It appears that seed values are taken from myelinated_axon_table, judging from the arguments in the extrude command.

· There are now notes in the new Geometric model on how to correct the indices in all_axon_table to match myelinated_axon_table.

SCIRun Data Analysis in Matlab

Set the anode and intensity values in analyze_output.m, then run it. This sets up the matrices, runs convert_sci_data.m and find_max_current_in_axon.m, and saves output to three text files. One important difference between the original method for compiling results and the new method used here is that the old method analyzed the entire mesh, while the new method only looks in the axoplasm. This is partly a side effect of inserting the cell membrane, but it also requires far less data to be moved around.

Old Instructions – probably ok to delete these

% load the variables from the scirun output file. If there is a read error in

% matlab, make sure it was exported in matlab format, and that you didn't get

% the scirun file on accident (both use the .mat extension)

load ../solution_data/050103/vector_ijk_anode1.mat

% give the variables better names

scirun_current_magnitude_anode1=i2;

clear i2

% repeat the above for the xyz location vectors (which is of size <#nodes>x3)

load ../solution_data/050103/vector_xyz_anode_xyz.mat

scirun_xyz=i1;

electrode_xyz=i3;

clear i1

clear i3

% Reshape the data from scirun to make a five-dimensional matrix

% first dimension is which tetvol in a given hexvol - five possibilities, so size=5

% second dimension is x-direction in original bitmap

% third dimension is y-direction in original bitmap

% fourth dimension is z-direction - size is equal to the number of extrusion layers

% fifth dimension is for direction-vector of current-density field. Three components

% to a 3d vector, so size=3

temp_linear_current=reshape(scirun_current_magnitude_anode1,[5,163,149,153,3]);

% use only one tetvol for each hexvol; this deletes 80% of the data

temp_linear_current=temp_linear_current(5,:,:,:,:);

% reshape the array to put it in four-dimensional form (otherwise, five dimensions

% of previous array will persist in matlab

scirun_current_magnitude_anode1=reshape(temp_linear_current,[163*149*153,3]);

% delete intermediate variable

clear temp_linear_current

% perform the same steps for the location-vectors that were done for the current-

% density vectors:

temp_linear_xyz=reshape(scirun_xyz,[5,163,149,153,3]);

temp_linear_xyz=temp_linear_xyz(5,:,:,:,:);

scirun_xyz=reshape(temp_linear_xyz,[163*149*153,3]);

% clear all temporary variables

clear temp*

% find the maximum current density in the z direction (parallel to the fiber axis)

% by using a custom function, the xyz location vector field, and the current density

% vector field

max_current_density_dot_z=convert_sci_data(scirun_xyz,scirun_current_magnitude_anode1);

% make sure you have all the right variables saved in your geometry data file

whos -file ../geometry_data/050203/nerve_model_050203.mat

% load the all_axoplasm_table and myelinated_axon_table from the relevant data

% files that were made in the model generation steps. These will allow us to convert

% data from the xyz location vectors that matlab gives us back to the xyz locations that

% correspond to our original bitmap.

load('../geometry_data/050203/nerve_model_050203.mat','all_axon_table','axoplasm_table');

load('../geometry_data/050203/nerve_model_050203.mat','myelinated_axon_table');

% find the maximum current in each axon in the all_axon_table by using the matrix of

% maximum current density data and a custom function.

max_current_in_axon=find_max_current_in_axon(all_axon_table,axoplasm_table,max_current_density_dot_z);

% save the resulting table, which has (among other things) the axon number, x and y locations, and maximum current

% density in the z direction

save('../data_analysis/050203/excel_max_current_in_axon_050203.txt','max_current_in_axon','-ASCII');

save('../data_analysis/050203/excel_all_axon_table_050203.txt','all_axon_table','-ASCII');

save('../data_analysis/050203/excel_myelinated_axon_table_050203.txt','myelinated_axon_table','-ASCII');

save ../solution_data/050203/analysis_050203.mat