#Workflow example#
import os
GUIinclude = True
from datetime import date
from carputils import settings
from carputils import tools
from carputils import mesh
from carputils import testing
from carputils import ep
from carputils.carpio import txt
from carputils.mesh import Block, generate
import numpy as np
from numpy import array as nplist
def parser():
parser = tools.standard_parser()
group = parser.add_argument_group('experiment specific options')
group.add_argument('--ionic model',
type = str,
default = 'courtemanche',
choices = ['tenTusscherPanfilov','HodgkinHuxley', 'DifrancescoNoble'],
help = 'pick ionic model, use bench to list options')
group.add_argument('--sourceModel',
type = str,
default = 'monodomain',
choices = ['monodomain','bidomain'],
help = 'Pick a domain model, mono- or bi- (default: monodomain)')
group.add_argument('--duration',
type = float,
default = 1000.,
help = 'Duration of simulation in [ms] (default: 1000.)')
group.add_argument('--stim-strength',
type = float,
default = 30.,
help = 'pick transmembrane current stimulus strength in [uA/cm^2] = [pA/pF] considering fixed Cm (default: 30.)')
group.add_argument('--stim-dur',
type = float,
default = 2.,
help = 'pick transmembrane current stimulus duration in [ms] (default: 2.)')
return parser
#Execution and file outputs
def jobID(args):
today = date.today()
return '{}_basic_{}'.format(today.isoformat(), args.duration)
@tools.carpexample(parser, jobID)
def run(args, job):
#Define the mesh#
x = 5.0
y = 5.0
z = 1.0
res = 0.1
#Defining the mesh, give it a centre, sizes in each dimension, a resolution and a mesh type
geom = mesh.Block(centre = (0.0,0.0,0.0), size = (x,y,z), resolution = res, etype = 'tetra')
#Define regions within the mesh
reg1 = mesh.BoxRegion((-x/2,-y/2,-z/2), (0, y/2, z/2), tag = 10)
geom.add_region(reg1)
reg2 = mesh.BoxRegion((0,-y/2,-z/2), (x, y/2, z/2), tag = 11)
geom.add_region(reg2)
#shifts origin to the bottom left corner
geom.corner_at_origin()
#sets fibres
geom.set_fibres(0,0,0,0)
meshname = mesh.generate(geom)
# query for element tags
_, etags,_ = txt.read(meshname + '.elem')
etags = np.unique(etags)
IntraTags = etags[etags != 0] # element labels for extracellular grid
ExtraTags = etags.copy() # element labels for intracellular grid
#Defining ionic models and conductivities for each region#
#ionic models
imp_reg = ['-num_imp_regions', 2,
'-imp_region[0].im', 'courtemanche',
'-imp_region[0].num_IDs',1,
'-imp_region[0].ID[0]', "10"
'-imp_region[1].im', "courtemanche",
'-imp_region[1].im_param', "g_Kr*2.,g_Ks*2.",
'-imp_region[1].num_IDs',1,
'-imp_region[1].ID[0]', "11"
#conductivities
g_reg = ['-num_gregions', 2,
'-gregion[0].num_IDs', 1,
'-gregion[0].ID[0]', 10,
'-gregion[0].g_il', 0.174,
'-gregion[0].g_it', 0.174,
'-gregion[0].g_in', 0.174,
'-gregion[0].g_el', 0.625,
'-gregion[0].g_et', 0.625,
'-gregion[0].g_en', 0.625,
'-gregion[1].num_IDs', 1,
'-gregion[1].ID[0]', 11,
'-gregion[1].g_il', 0.074,
'-gregion[1].g_it', 0.074,
'-gregion[1].g_in', 0.074,
'-gregion[1].g_el', 0.025,
'-gregion[1].g_et', 0.025,
'-gregion[1].g_en', 0.025]
#Defining the stimulus#
#Electrode definition
E1_lower_bound = nplist([ 0, 0, z])
E1_upper_bound = nplist([ 2*res, 2*res, z])
electrode = mesh.block_region(geom, 'stimulus', 0, E1_lower_bound, E1_upper_bound, False)
electrode = electrode
stim = ['num_stim', 1,
'-stimulus[0].stimtype', 0,
'-stimulus[0].strength', 2500.0,
'-stimulus[0].duration', 2.0,
'-stimulus[0].start', 0,
'-stimulus[0].npls', 1]
#Define simulator outcomes#
# Gives model type options
Src = ep.model_type_opts(args.sourceModel)
#dt gives time resolution for calculations, spacedt gives the interval for which values are written out timedt givs interval for the command line to give updates
num_par = ['-dt', 100]
I0_par = ['spacedt', 0.1,
'timedt', 5.0]
#The command line inputs
cmd = tools.carp_cmd()
cmd += imp_reg
cmd += g_reg
cmd += stim + electrode
cmd += num_par
cmd += I0_par
cmd += Src
simID = job.ID
cmd += ['meshname', meshname,
'-tend', args.duration,
'-simID', simID]
if args.visualize:
cmd += ['gridout_i', 4,
'gridout_e', 4]
job.carp(cmd)
if __name__ == '__main__':
run()
