This example illustrates the use of the Noble 98 cellular electrophysiological
model in a 3-dimensional mesh. The geometry used is a two element wedge taken
from the LV free wall of the rectangular-Cartesian porcine heart model. The
mesh looks like this (with fibres)...
The model is solved using implicit Adams-Moulton integration with grid-based
FEM. The tissue is stimulated at the grid points along the lower endocardial
surface.
# Example b211v
# Noble 98 activation in a 3d LV free wall wedge.
{
eval { require Time::HiRes };
if( $@ )
{
print "warning: $@";
}
else
{
import Time::HiRes qw(time);
}
}
sub cpu_time()
{
# Add user/system parent/child cpu_times.
my $result = 0;
map { $result += $_ } times;
$result;
}
sub print_timing($$$)
{
my ($message,$wall,$cpu) = @_;
printf "%s: %.3fs wall, %.3fs cpu.\n", $message, $wall, $cpu;
}
sub print_memory_use()
{
# Different operating systems have different data and ps uses different
# arguments so try a series of commands in an attempt to collect the
# information.
system qw(ps -o vsz -p), $$;
system qw(ps -l -p), $$;
system qw(ps v), $$;
}
$output = ".";
$format = "binary";
# If the example path is not set, default to current directory
if (! defined $example) {
$example = "./";
}
# Drop off the trailing / in the example path
$chopped = chop $example;
if ($chopped ne "/") {
$example .= $chopped;
}
set echo;
# Define the model size
fem define parameters;r;$example/minimal;
# Define the geometry of the model
fem define coordinate;r;$example/emech;
fem define node;r;$example/wedge;
fem define base;r;$example/emech-cubic-linear;
fem define element;r;$example/wedge;
fem group element 10,29 as ALL_ELEMENTS;
fem define fibre;r;$example/wedge;
fem define element;r;$example/wedge fibre;
# Define the grid scheme to use
fem define grid;r;$example/cell-60x60x40;
fem update grid geometry;
fem update grid metric;
fem group grid external as boundary;
# Create the groups to use as the stimulus site
fem group grid line 147482 as group1 xidirn 1 positive;
fem group grid line 147423 as group2 xidirn 1 positive;
fem group grid line 147364 as group3 xidirn 1 positive;
fem group grid line 151022 as group4 xidirn 1 positive;
fem group grid line 150963 as group5 xidirn 1 positive;
fem group grid line 150904 as group6 xidirn 1 positive;
fem group grid line 154562 as group7 xidirn 1 positive;
fem group grid line 154503 as group8 xidirn 1 positive;
fem group grid line 154444 as group9 xidirn 1 positive;
fem group grid line 158102 as group10 xidirn 1 positive;
fem group grid line 158043 as group11 xidirn 1 positive;
fem group grid line 157984 as group12 xidirn 1 positive;
fem group grid grid group1,group2,group3,group4,group5,group6,group7,group8,group9,group10,group11,group12 as stimulus_site;
# Define the Noble 98 electrical cellular based model, using
# grid-based FEM
fem define equation;r;$example/n98;
# Define the cellular material parameters
fem define cell;r;$example/n98;
# and their spatial variation
fem define material;r;$example/n98 cell;
# Define the distributed material parameters
fem define material;r;$example/n98;
fem update grid material;
# Define boundary conditions
fem define initial;r;$example/n98;
# Define an implicit integration scheme using the
# adaptive Adams-Moulton ODE integrator
fem define solve;r;$example/implicit_grid_based_fem;
# Export the geometry and grid point numbers
fem export node;$output/wedge as wedge;
fem export element;$output/wedge as wedge;
fem export element;$output/numbers as wedge grid_numbers;
($wall_setup,$cpu_setup) = (time, cpu_time);
print_timing 'setup', $wall_setup - $^T, $cpu_setup;
# Initialise the solution
fem solve to 0;
($wall_matrix,$cpu_matrix) = (time, cpu_time);
print_timing 'matrix generation', $wall_matrix - $wall_setup, $cpu_matrix - $cpu_setup;
# Get the index of the membrane potential variable
fem inquire cell_variable Vm return_variables VM_ARRAY,VM_IDX;
# and create a history file to store its value through time
fem open history;$output/wedge write variables $VM_ARRAY niqslist $VM_IDX $format;
($wall_output,$cpu_output) = (time, cpu_time);
print_timing 'output open', $wall_output - $wall_matrix, $cpu_output - $cpu_matrix;
# The time loop
for $time ( 0..2 )
{
# Solve the model
fem solve restart to $time;
($wall_solve,$cpu_solve) = (time, cpu_time);
print_timing "solve $time", $wall_solve - $wall_output, $cpu_solve - $cpu_output;
# Write the solution to the history file
fem write history time $time variables $VM_ARRAY $format;
# and export the potential field
fem export element;$output/potential_$time field as wedge;
($wall_output,$cpu_output) = (time, cpu_time);
print_timing "output $time", $wall_output - $wall_solve, $cpu_output - $cpu_solve;
}
fem close history $format;
print_memory_use;
# Generate a CMISS signal file
#fem evaluate electrode;$output/Vm history yqs iy $VM_IDX $output/wedge from grid $format;
# and export it to UnEMAP
#fem define export;r;$example/cell;
#fem export signal;$output/Vm electrode signal $output/Vm;
