This document briefly explains the modeling process, from microscope slide to experimental results
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.
=====================================================
% 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
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.
Once the mesh is in SCIRun format there are several steps to complete the biophysical model.
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
· 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.
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.
% 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