Abstract
As part of a 3wk intersession workshop funded by a National Science Foundation Expeditions in Computing award, 15 undergraduate students from the City University of New York^{1} collaborated on a study aimed at characterizing the voltage dynamics and arrhythmogenic behavior of cardiac cells for a broad range of physiologically relevant conditions using an in silico model. The primary goal of the workshop was to cultivate student interest in computational modeling and analysis of complex systems by introducing them through lectures and laboratory activities to current research in cardiac modeling and by engaging them in a handson research experience. The success of the workshop lay in the exposure of the students to active researchers and experts in their fields, the use of handson activities to communicate important concepts, active engagement of the students in research, and explanations of the significance of results as the students generated them. The workshop content addressed how spiral waves of electrical activity are initiated in the heart and how different parameter values affect the dynamics of these reentrant waves. Spiral waves are clinically associated with tachycardia, when the waves remain stable, and with fibrillation, when the waves exhibit breakup. All in silico experiments were conducted by simulating a mathematical model of cardiac cells on graphics processing units instead of the standard central processing units of desktop computers. This approach decreased the run time for each simulation to almost real time, thereby allowing the students to quickly analyze and characterize the simulated arrhythmias. Results from these simulations, as well as some of the background and methodology taught during the workshop, is presented in this article along with the programming code and the explanations of simulation results in an effort to allow other teachers and students to perform their own demonstrations, simulations, and studies.
 ventricular fibrillation
 atrial fibrillation
 graphics processing unit simulations
the 2011 atrial fibrillation workshop was one of a series of annual workshops affiliated with the Computational Modeling and Analysis of Complex Systems (CMACS), a multiinstitutional (8) and multiprincipal investigator (19) project led by Edmund Clarke. CMACS is funded by a National Science Foundation Expeditions in Computing award. The objective of the workshops is to develop the scientific interest and skills of students from urban minorityserving institutions and especially to motivate them to study the kinds of computational modeling techniques and applications used and developed in the project. Each year, the workshop is held in January at Lehman College in the Bronx, New York. Students from all parts of The City University of New York (CUNY), in which Lehman College is only one of many colleges, are invited to apply for the workshop; applicants are admitted based on grades and recommendations. The first workshop was held in January 2010 and addressed the role of cellular signaling pathways in the development of pancreatic cancer. Instruction at the 2010 workshop was provided by faculty members from Lehman College (N. Griffeth), the University of Pittsburgh (James Faeder), Carnegie Mellon (Christopher Langmead), and New York University (Bud Mishra) together with several CUNY graduate students (Loes Olde Loohuis, Ziping Liu, and Fred Dieckamp). For the 2011 workshop, the instructors were from Cornell University (F. H. Fenton and R. F. Gilmour, Jr.), Stony Brook (E. Bartocci, S. A. Smolka, and J. Glimm), and Lehman College (N. Griffeth) together with several Cornell and CUNY graduate students (R. Singh, F. B. von Stein, T. GrossoApplewhite, K. Zhao, J. Rogers, and A. Wolinetz).
The Workshop
The Atrial Fibrillation workshop consisted of three related weeklong units, each with a series of lectures and laboratory exercises, conducted 5 days/wk for 5 h/day. The teaching philosophy of the workshops is to engage the students in discussion, group work, and activities more than lectures. The workshop was organized around alternating lectures with laboratories and group activities to implement this philosophy.
In the first week, F. H. Fenton and R. Singh gave a series of lectures and laboratory exercises on excitable systems in biology (41, 42, 66, 72), physics (39, 45, 55), chemistry (6, 8, 78), and mathematics (26, 43, 48). This background provided students with a robust understanding of excitable systems and their relationship to complex systems and chaos (35, 71, 77). Other material in the first week emphasized the role of both mathematical and computational frameworks in the study of cardiac arrhythmias (16, 31). The motivation for computational and mathematical approaches in the study of cardiac arrhtyhmias was further conveyed to the students by explaining how the field of cybernetics, founded by Norbert Wiener, was a result of his studies not only on antiaircraft defense but also of cardiac arrhythmias such as fibrillation (75) at the beginning of the 1940s. Emphasis was given to inquirybased teaching. As a result, for the first part of the week, students experimented with density (saline) oscillators (51) and chemical oscillators (6, 8, 52, 78) to relate the nature of the oscillations in each system with those of cardiac action potentials (APs) and to quantify in a simple way the similarities between the complex dynamics in these systems and those of cardiac arrhythmias. The second part of the week concentrated more on how the membrane potential of cells such as neurons and cardiomyocytes can be represented mathematically. Various models (26) with different levels of biophysical detail were simulated and studied using interactive Java applets (24) (created by E. M. Cherry, A. B. Filipski, and F. H. Fenton). At the end of the first week, R. F. Gilmour, Jr., gave a lecture on experimental methods that illustrated how the electrophysiological properties F. H. Fenton had introduced in the context of computational models could be observed and measured experimentally.
The second week began with an exercise (designed by T. GrossoApplewhite) with two goals in mind. The first, overt goal was to provoke the students to think critically about the value of models and what characteristics make them valuable. The second, unstated goal was motivated by observations from the first workshop that some of the student groups tended to degenerate into three people working side by side instead of collaboratively. This goal was to acquaint the students with each other and to push them out of their usual mode of individual student work, encouraging them instead to contribute actively to the group's progress. The exercise asked students to first list the characteristics of good models, individually, and then to prioritize them as a group. Subsequent group work validated the success of this exercise.
N. Griffeth and K. Zhao lectured on various numerical methods used to analyze the mathematical models presented during the first week, and J. Rogers presented a tutorial on Compute Unified Device Architecture (CUDA) programming (62). CUDA provides a programming interface for the highly parallel graphics chips (graphics processing units) from NVIDIA, which are now standard in many laptops and desktops. Because CUDA programming uses the multiprocessing capabilities of the NVIDIA chips (3, 62), it provides much faster calculations than those obtained using the central processing units. For the workshop, J. Rogers and A. Wolinetz (both Master's degree students in the Lehman College Computer Science program) made available two computers that they had built around these highperformance NVIDIA graphics chips and supported their use by the undergraduate students. At the end of the second week, E. Bartocci (Stony Brook University) gave an indepth series of lectures on the implementation of an electrophysiological model of a cardiac cell (9) in a CUDA program.
During the final week, E. Bartocci and the students modified the E. Bartocci, E. M. Cherry, and F. H. Fenton CUDA program to simulate twodimensional idealized sections of cardiac tissue and to examine a set of five different parameter spaces in detail to identify how certain changes in physiology affect the stability and dynamics of reentrant waves (spiral waves). The students subsequently prepared presentations with their results and interpretations, which were given on the final day of the workshop. All lectures and exercises can be found at http://www.lehman.edu/academics/cmacs/outline.php.
The Saline Oscillator
As a first approach to introducing students to the concepts of nonlinear dynamics and chaos necessary to understand cardiac dynamics (35, 71), in general, and arrhythmias, in particular, from both a physical and mathematical point of view, F. H. Fenton and R. Singh prepared a series of exercises using several setups of saline oscillators. A density (in this case, saline) oscillator is a very simple and inexpensive, yet extremely interesting, system discovered by accident in the late 1960s by Martin Seelye (51), where a container with a higherdensity fluid is immersed into another fluid with a lower density with a narrow orifice connecting the two containers (see Fig. 1A). The difference in densities, combined with gravity, leads to a RayleighTaylor instability (37, 64) that results in a rhythmic change in flow of fluid from one container to the other that can last for many hours (12, 39, 58). Although the phase transition dynamics are complicated (59), qualitatively the denser fluid has a tendency to fall and the less dense fluid has a tendency to rise. Separating the two solutions via hydrostatic pressure alone is a difficult balancing act. Typically, a small amount of solution will escape in one direction or the other, creating a jet. Over some period of time, this jet eventually will taper off as too much liquid accumulates on the receiving container; the imbalance in hydrostatic pressure then drives a jet in the other direction. The phenomenon continues as long as the density difference across the orifice is large enough to sustain it. This flow or leak of saline water (when going from the upper chamber to the lower chamber) and fresh water (when going from the lower chamber to the upper chamber) is clearly visible (see the movie at http://TheVirtualHeart.org/CMACS/salineoscillator.mov) as a laminar flow through the orifice and as a buoyant jet at the orifice exits. Furthermore, the flux across the orifice can be easily quantified using two electrodes and an oscilloscope. By placing one electrode in line with the orifice and the other far away from it, oscillatory changes in the boundary layer of the electrode in line with the orifice can be used to measure a change in voltage of ∼100 mV that quantifies the presence and absence of a jet (73). This oscillating voltage signal is not only qualitatively similar to an AP from an autooscillatory sinoatrial (SA) node cell but, in some respects, is also qualitatively similar to cardiomyocyte APs in general (see Fig. 1B). In both cardiac cells and the oscillator, it is possible to provide an initiating stimulus [a small injection of fresh water in the outer container for the oscillator (39)], and, by constant periodic stimulation, it is possible to induce complex rhythms such as alternans (see Fig. 1C) (40, 60, 74) and Wenckebach rhythms (23, 34, 69). Therefore, the oscillator is an excellent system to provide students with a handson experience demonstrating how excitable systems behave and how the oscillator, in particular, shares many dynamic similarities with cardiac cells.
For the experiments, students were divided into four groups, each of which studied the effects of changing four physical parameters on the system's dynamics. In particular, each team tried several different orifice diameters, orifice lengths, density differences, and chamber sizes and recorded the changes in the oscillation period. After cycling through each of the four exercises, each team was asked to plot their results, compare their findings with those of the other groups, and supply an explanation for the observed trends based on their physical intuition. Afterward, a detailed explanation of the results was given in a lecture. Figure 1D shows one example in which the orifice length was changed. For the orifice length chosen in our setup, the time scales of the laminar flow were larger than the time scales corresponding to the transitions in the flow reversal; as a result, the flow rate through the orifice was assumed to be largely governed by the HagenPoiseuille law (44). An increase in orifice length led to a linear decrease in flow rate and, in turn, to an approximately linear increase in the period of oscillation.
Chemical Oscillators
Although the saline oscillator allowed the students to understand some of the concepts of complex systems that are also found in cardiac dynamics, such as oscillations between unexcited and excited states (polarizationdepolarization) along with alternans and higherorder periodicities (Wenckebach rhythms), the dynamics were confined to a single cell or oscillator, so it was not possible to study spatial effects. Chemical oscillators, however, can be used to easily demonstrate some of the additional levels of complexity that arise when coupling in space is considered. The BriggsRauscher (BR) reaction (8) is an oscillator similar to the saline oscillator in that the system continually oscillates, in this case, between two colors (see Fig. 2). Just as in the case of the SA node and the saline oscillator, the period in this system also can be changed by modifying some of the chemical constituents (11, 33), in this case, the amount of hydrogen peroxide used in the recipe. As long as the liquid remains homogenous (i.e., continually being swirled), the change in color occurs everywhere at the same time. However, if the liquid is not perturbed in this manner, small heterogeneities in the mixture will result in some regions oscillating out of phase, so that some locations will change color before others. This is the principle behind propagating waves in space, where only part of the domain is excited, and the excited state can advance by exciting the unexcited (quiescent) part of the domain.
To illustrate the propagation of waves, F. H. Fenton and F. B. von Stein created a set of experiments for the students using the BeluosovZhabotinsky (BZ) reaction (6, 52, 78), which is a system that oscillates between red and blue at a much slower time scale than that of the BR reaction. In effect, the initiation and visualization of wave propagation (22, 63), as shown in Fig. 3A, is much easier. Each student was given a petri dish and instructed to pour in enough BZ solution using a pipette to form a thin (quasitwodimensional) layer 1 mm deep (Fig. 3B). Students then observed the initiation of target waves due to small localized changes in concentrations. The sensitivity to initial conditions (butterfly effect) (36, 67) of chaotic systems was then explained to the students as they observed that each petri dish produced a completely different pattern even though all petri dishes were started in a similar way. They understood that even very small differences in initial conditions could make the dynamics of two instantiations diverge and produce completely different results (patterns in this case) over time. Students were encouraged to mix the petri dish many times to restart the experiment and observe the different patterns created every time.
To mimic the propagation of electrical waves in the heart, which are normally produced by the SA node or ectopic foci, students used silver wire to touch the BZ liquid (exciting the system), thereby generating propagating waves. The silver metal in the oxidizing solution forms Ag^{+}, which reacts with Br^{−} to form the precipitate AgBr. This causes an excitation since Br^{−} is an inhibitor of autocatalysis (68). This variation in concentration starts to propagate by diffusion, giving rise to a propagating wave, where the size of the area of the propagating wave is a function of the diameter of the silver wire and the time that it is in contact with the solution (38).
Students then proceeded to generate spiral waves using a method similar to the initiation protocols for cardiac tissue known as the S1S2 protocol (32) or pinwheel protocol proposed by Winfree (77) modified for the BZ media (30, 38). This method generates a first wave by an initial stimulus S1 (in this case, by touching the silver wire to one section of the petri dish), after which a second stimulus S2 (again using the silver wire) is induced somewhere behind the first wave. If this second wave is at just the right distance from the back of S1, it will find a vulnerable window (65, 70) where part of the system is excitable and part is refractory, so that the second wave can propagate in only one direction, initially. As the medium immediately behind the first wave recovers, the new wave front curls, thereby initiating two counterrotating spiral waves (see Fig. 3, C and D). The generation of these spiral waves is equivalent to an induction of tachycardia in cardiac tissue, because the spiral waves have a rotation period faster than that of the natural oscillation, thereby allowing them to take control of the system at this faster period. In cardiac tissue, the natural oscillation corresponds to the activations produced by the SA node (sinus rhythm) and the faster period of tachycardia is produced by a single spiral wave or one or two pairs of counterrotating spiral waves.
The Student Research Experience
The final step of the workshop was to engage the students in a handson research experience. In this section, we explain the approach to engaging the students, including the elements that made this part of the workshop successful, as well as the accessibility of Java applets to develop intuition about the processes being studied and ongoing discussions about the significance of the work.
Research objective.
The objective of the studentinvolved research was to study the dynamics of a cardiac cell model in different parameter regimes. Students were instructed that spiral waves in cardiac tissue can exhibit different dynamics depending on the physiology and structure of the tissue (2, 13, 18, 20, 25, 46, 47, 76), similar to the earlier laboratory exercise with the BZ reaction (5, 7, 49, 54). Spiral waves can rotate around small obstacles or have a very small radius of rotation (see Fig. 4A), or they can have a long linear core trajectory (see Fig. 4B). In many cases, initiated spiral waves are unstable and quickly break up into multiple spiral waves (see Fig. 4C). In between these limits, there is a wide range of other possible trajectories. For chemical reactions such as the BZ reaction and for different mathematical models of cardiac cells, there have been a number of studies aimed at characterizing spiral wave dynamics as a function of model parameters (2, 13, 18, 20, 25, 46, 47, 76) to establish which physiological changes are more or less proarrhythmic.
Student preparation.
During the first 2 wk, the lectures on excitable systems by F. H. Fenton and the laboratory experiments with saline oscillators and the BZ reaction gave the students a basic understanding of how the dynamics of complex systems can be modeled mathematically. As part of the early exercises, students were given time to work with a series of interactive Java applets (24) to explore the dynamics of various cardiac cell models in a single cell and in tissue (http://TheVirtualHeart.org/applets.html) and were encouraged to use the interactive Java applet for the simulation model in a single cell. Tutorials on numerical integration by K. Zhao and CUDA programming by J. Rogers prepared the students for the implementation of the cardiac cell model in CUDA by E. Bartocci, E. M. Cherry, and F. M. Fenton.
Group composition.
In the third week, students were divided into five groups. Each group contained three students and included at least one student with courses in biology and one student with courses in computer science. Almost all of the students had taken calculus, so that was not a criterion. The objective for each group was to characterize the dynamics of spiral waves as a function of different physiological parameters to understand the effects of physiology on the simulated arrhythmias.
Methods.
Students used the CUDA program downloadable from http://TheVirtualHeart.org/CMACS/4vmCUDA.tar to simulate the dynamics for many different values of the parameters. The in silico model used for this study is a reduced ionic model that reproduced the AP of a human epicardial cell (9), including experimentally measured AP shape, upstroke amplitude, and rate of rise, threshold of excitation, adaptation to changes in cycle length (i.e., the restitution curve), and conduction velocity. Briefly, in this model, the cardiac AP is modeled by three main currents: a fast inward current responsible for depolarization (positive change in voltage) that can be associated with the fast Na^{+} current, a slow outward current that is responsible for repolarization (negative change of voltage) and is analogous to a timedependent K^{+} current, and a slow inward current analogous to a total Ca^{2+} current that balances the slow outward current after depolarization to produce the plateau of the AP. Some of the model parameters can be interpreted as analogs of components of other cardiac electrophysiology models of comparable complexity, such as the BeelerReuter (4) and LuoRudy I models (50). An interactive Java applet for this model in a single cell that allows variation of any of the 28 parameters of the model to analyze how the AP characteristics change as a function of the different parameter values, including those studied by the students and presented below, can be run from http://TheVirtualHeart.org/java/4vn/fourv0d.html.
In each group, students initiated a spiral wave in the computational model in a manner similar to that used in the BZ reaction experiments (i.e., by an S1S2 premature stimulation) and changed two specific sets of parameters to characterize the resulting dynamics and construct a parameter space plot. The parameters and their values are shown in Table 1. A total of 212 simulations were performed, but because the code ran essentially in real time (3), students were able to observe in a fast, easy, and interactive way how the dynamics of the reentrant waves varied as a result of parameter changes. The CUDA code is implemented with a tipfinding algorithm (27) that plots the trajectory of the spiral wave tips along with the spiral wave as the spiral evolves in time. The numerical integration used conforms with the benchmark of cardiac cell modeling (56); however, as the parameters are changed and the upstroke becomes steeper or the AP becomes very short, larger temporal or spatial resolution is required, respectively. The temporal resolution was sufficient for all the cases studied. In most cases, we used a standard fivepoint discretization of the Laplacian. However, for values of the time constant of the slow outward current during repolarization (τ_{so1}) 10, this combination of discretization scheme and spatial resolution is insufficient and can lead to pinning to the lattice and the appearance of corners within the spiral waves (16). To ensure consistency in our comparisons, we opted to keep the same domain size fixed for all simulations and for τ_{so1} < 10 instead used a ninepoint Laplacian discretization whose leading error term is rotationally symmetric and thus minimizes lattice artifacts (16, 56).
Student Results
Students spent 5 days preparing and running simulations. Intermittently, F. H. Fenton engaged in discussion with the students to explain the results as they developed. Explanations are included below with each group's results. We believe that this kind of interaction with the students emulates the progress of a professional research project and is important to maintaining their enthusiasm for the work.
Group 1 varied the following parameters: the time constant to inactivate the slow inward Ca^{2+} current (τ_{w}^{+}) and the time constant to recover from the inactivation of the slow inward Ca^{2+} current (τ_{w}^{−}). Both parameters affect the dynamics of the inactivation gate variable w, which for the model used (9) is analogous to the voltagedependent inactivation gate variable (f) in the BeelerReuter and LuoRudy I models. However, for this model, the time constant for inactivation is a step function, with τ_{w}^{+} active during depolarized potentials and τ_{w}^{−} active during polarized potentials, which can be interpreted as the time constant to inactivate τ_{w}^{+} and the time constant to recover from inactivation τ_{w}^{−} of the Ca^{2+} current. As the students could see from the Java applet, an increase (decrease) in τ_{w}^{+} effectively increases (decreases) the AP duration (APD). The parameter τ_{w}^{−}, on the other hand, dictates the time it takes after an activation for a subsequent APD to reach the maximum APD, thereby affecting restitution properties (10, 14, 19, 53, 57, 61), where the larger the time constant, the longer it will take to reach the maximum APD.
Figure 5A shows the spiral wave profiles for the 64 different sets of parameters, and Fig. 5B shows their corresponding tip trajectories. Note that for better visualization of the trajectories, they all have been normalized to fit the same area; the real size can be estimated by observing the spiral waves. For small values of τ_{w}^{+}, the fast inactivation of the slow inward Ca^{2+} current results, as expected, in short wavelengths, because the Ca^{2+} current is decreased, thereby reducing the plateau and leading to short APDs. As τ_{w}^{+} increases, inactivation takes longer, resulting in increased Ca^{2+} current that produces longer APDs and longer wavelengths that yield thicker spiral waves. A clear transition can be seen between τ_{w}^{+} = 150 ms and τ_{w}^{+} = 200 ms; in this transition, the spiral wave trajectory changes from hypocycloidal to linear.
As τ_{w}^{−} increases, the time to recover from inactivation increases, which results in a large dispersion when the APD becomes longer (i.e., large values of τ_{w}^{+}). Increased dispersion is known to be proarrhythmic (1, 17, 60, 74), in this case, by alternans (a rhythm previously shown to the students in the saline oscillator). In fact, as shown in Fig. 5A, breakup of spiral waves is obtained for larger values of both time constants. Figure 5B shows this region, where no lasting single spiral wave trajectory can be drawn, in dark gray, whereas the lighter gray colors indicate the transition region between stable and unstable spiral waves where quasibreakup (transient breakup of a single spiral wave) occurs.
Group 2 varied the following parameters: the time constant for the slow inward Ca^{2+} current (τ_{si}) and the time constant for the Na^{+} or fast inward current (τ_{fi}). The parameter τ_{si} is the inverse of the maximum Ca^{2+} conductance or slow inward current, analogous to 1/g_{si} in the BeelerReuter and LuoRudy I models, and τ_{fi} is the inverse of the maximum conductance for the Na^{+} or fast inward current, analogous to 1/g_{Na} in the BeelerReuter and LuoRudy I models. The larger (smaller) τ_{si}, the smaller (larger) the Ca^{2+} current and the shorter (longer) the APD; the larger (smaller) τ_{fi}, the lower (higher) the excitability of the tissue, as the amount of Na^{+} current entering the cell during depolarization decreases (increases). Again, students could verify these behaviors using the Java applet.
Figure 6 shows that the larger the APD or wavelength (τ_{si} < 2.0), the longer the spiral wave tip needs to travel before finding excitable tissue to make a turn, resulting in long linear cores (21, 47). As excitability decreases, the tip takes longer to make a turn, resulting in more precession in the linear core trajectory. As τ_{si} increases, the wavelength decreases, and the spiral wave meanders more for low excitability, but as the excitability increases (smaller τ_{fi}), the faster turns make the tip of the spiral collide with its own wave back, producing quasibreakup. Figure 6B shows the different levels of quasibreakup and meander as a function of the wavelength (determined here by τ_{si}) and excitability (determined here by τ_{fi}) of the simulated tissue.
Group 3 varied the following parameters: τ_{so1} and τ_{si} (described above). The parameter τ_{so1} is the inverse of the maximum conductance of the slow outward current during depolarization, so it is similar to 1/ḡ_{x1} in the BeelerReuter model and 1/Ḡ_{K}X_{i} in the LuoRudy I model; τ_{si} is as described above.
The larger (smaller) τ_{so1}, the smaller (larger) the effective K^{+} current and the larger (smaller) the APD. For small values of τ_{si}, the repolarization current is so strong that the APD becomes as small as that of a mouse or rat, resulting in very thin spiral waves similar to those obtained with the BZ reaction, as shown in Fig. 7. As the repolarization current decreases (larger τ_{so1}), the AP plateau duration increases, and the wavelength of the spiral increases. As the Ca^{2+} current increases (smaller values of τ_{si}), the tip meanders less, and the linear cores precess less. Figure 7 shows that for these sets of parameter values, intermediate values of the repolarization current (τ_{so1}) produce a quasibreakup regime. This occurs in the intermediate regime in between meander and linear cores, where complex hypermeandering trajectories are present. When hypermeander is complex, in some cases, the spiral wave tip can turn so fast that a Doppler shift effectively increases the frequency of the spiral wave locally to the point that it reaches the refractory period, producing temporary (nonsustained) breakup that subsequently heals.
Group 4 varied the following parameters: τ_{so1} and τ_{fi} (both described previously). Figure 8 shows that as the excitability decreases and the repolarization current is quite strong, no spiral waves can be induced because there is not enough electrotonic current to produce an excitation to sustain wave propagation. As in the previous case, a large repolarization current results in spiral waves with small wavelengths; however, as the excitability increases (smaller τ_{fi}), the tip of the spiral wave can turn faster, as it needs to excite fewer neighbors to propagate (28, 47), thereby resulting in a faster rotation frequency and tighter spiral waves compared with those shown in Fig. 7. However, as with the results from group 3, there are values of depolarization and hypermeandering trajectories for which transient breakup by the Doppler effect occurs.
Group 5 varied the following parameters: the time constant from inactivation of the slow inward Ca^{2+} current (first part) (τ_{w1}^{−}) and the time constant of slow inward Ca^{2+} current (second part) (τ_{si2}). As in group 1, τ_{w1}^{−} represents the recovery from inactivation for the slow inward current, and, as in groups 2 and 3, τ_{si} is the time constant for the slow inward Ca^{2+} current. In the minimal model (9), these parameters can be piecewise functions between two values (e.g., τ_{w1}^{−} and τ_{w2}^{−}, and τ_{si1} and τ_{si2}, respectively); in this case, only τ_{w1}^{−} and τ_{si2} were varied while keeping the other two values constant at the values given in Ref. 9. τ_{w1}^{−} can produce large regions of dispersion when the APD is long, as in the case of small values of τ_{si2}. Because no parameter related to the upstroke (fast inward current) of the model has been changed, the excitability remains relatively high, resulting in relatively thick spiral waves (see Fig. 9A) that exhibit fast turns. As shown in Fig. 9B, there is a large region of spiral wave breakup and quasibreakup (shown in grayscale). In this case, the mechanisms described in Fig. 5 and those in Figs. 6⇑–8 are mixed. Most of the large dispersion and alternans breakup occurs within the region delimited by the solid line, whereas the Doppler effect breakup occurs within the region delimited by the dashed line.
Conclusions
The goals of the workshop were to cultivate student interest in the area of computational modeling and analysis of complex systems and to provide students early in their academic careers with a handson research experience. This article is itself evidence of their engagement in a handson research experience, and student feedback, as shown in Table 2, indicated that student interest was captured effectively.
Using simpler systems such as the saline oscillator and chemical oscillators, the students not only gained firsthand experience in experimental studies but also worked individually and as teams to answer questions that helped them understand the dynamics of complex systems and to participate actively in a research study that resulted in this report on which they appeared as coauthors. The experimental component therefore prepared them to develop a stronger understanding of how cardiac cells work, to learn how mathematical models are used for the study of arrhythmias, and to study how changes in physiological properties can alter the electrical dynamics of cardiac cells in tissue. From the five independent group studies, the students have shown that spiral wave dynamics and trajectories depend on model parameters. As has been observed in other models (2, 13, 18, 20, 25, 46, 47, 76), as the wavelength increases, two main effects can be seen: the likelihood of linear core trajectories increases and breakup becomes increasingly probable. Conversely, as the wavelength deceases, the likelihood of breakup decreases; however, the range of possible tip trajectories becomes much broader. In one particularly novel finding, students were able to demonstrate for the first time that, although previous studies of spiral waves with linear cores have always shown a precession that makes the angle of rotation vary, it is possible to obtain spiral waves with linear cores that do not precess for particular regions in parameter space with certain pairs of values for τ_{w}^{+} and τ_{w}^{−} (see Fig. 5B) or τ_{fi} and τ_{si} (see Fig. 6B).
We anticipate that some of the students from the workshop will continue working with the researchers leading the workshop over the summer, a result that also occurred in the last workshop, where two students obtained summer research jobs at Carnegie Mellon University.
GRANTS
This work was supported by National Science Foundation Grants CCF0926190 and CCF1018459 and in part by Air Force Office of Scientific Research Grant FA05500910481.
DISCLOSURES
No conflicts of interest, financial or otherwise, are declared by the author(s).
ACKNOWLEDGMENTS
The authors thank Humberto Arce and Hortencia Gonzales for the helpful advice in preparing the saline oscillator setups and Kenneth Showalter for comments on the chemical section. The authors also thank the Information Technology Department of Lehman College for access to the computer laboratory where most of the workshop took place and Brian Murphy for the use of CUDAready computers.
APPENDIX
The BR Reaction
A 50ml volume containing sodium iodate (328.5 mM) and sulfuric acid (130 mM) in water was mixed with an equal volume of a solution containing malonic acid (200 mM), manganese (II) sulfate monohydrate (26 mM), and soluble starch (2.5 mM) in water. A 100ml volume of 3% hydrogen peroxide was then added to create the reaction.
The BZ Reaction
A 50ml volume containing sodium bromate (250 mM) in water was mixed with an equal volume of malonic acid (310 mM) and sodium bromide (58 mM) in water. After the solution cleared, a 50ml volume of cerium (IV) ammonium nitrate (20 mM), sulfuric acid (2.7 M), and 3 ml ferroin (25 mM) was added to create the reaction.
The BZ Reaction for Waves
Seven milliliters of sodium bromate (0.5 M) in sulfuric acid (0.6 M) were mixed with 3.5 ml of malonic acid (0.5 M) in water and 1 ml of sodium bromide (0.97 M) in water. After the solution cleared, 1 ml of ferroin (25 mM) was added, and the solution was again mixed. A volume of this solution was transferred to a clean petri dish to form a thin layer 1–2 mm deep to observe spontaneous spiral wave formation.
Neutralization of the Chemical Reactions
BR reaction.
Neutralize the reduced iodide with the addition of 5 g sodium thiosulfate per liter of demonstration solution. Flush down the drain with water.
BZ reaction.
Neutralize with sodium bicarbonate. Flush down the drain with plenty of water.
CUDA Program
A copy of the CUDA code used by the students in this article can be downloaded from http://TheVirtualHeart.org/CMACS/4vmCUDA. To run it, one needs a computer with an NVIDIA card. You can download the CUDA drivers and toolkit from the NVIDIA web site: http://developer.nvidia.com/cudatoolkit40.
Install the CUDA drivers.
Install the CUDA toolkit.
Make ./cuda2DMinModel.
Below is a shortcut of the keyboard functions that can be used during simulations:
“P” key: pause the simulation
“R” key: restart the simulation
“Q” key: quit the simulation
“U” key: show the u field
“V” key: show the v field
“W” key: show the w field
“S” key: show the s field
“T” key: show or hide the tip trajectory
Footnotes

↵1 Brooklyn College, Hunter College, Lehman College, and Queens College.
 Copyright © 2011 The American Physiological Society
REFERENCES
 1.↵
 2.↵
 3.↵
 4.↵
 5.↵
 6.↵
 7.↵
 8.↵
 9.↵
 10.↵
 11.↵
 12.↵
 13.↵
 14.↵
 15.↵
 16.↵
 17.↵
 18.↵
 19.↵
 20.↵
 21.↵
 22.↵
 23.↵
 24.↵
 25.↵
 26.↵
 27.↵
 28.↵
 30.↵
 31.↵
 32.↵
 33.↵
 34.↵
 35.↵
 36.↵
 37.↵
 38.↵
 39.↵
 40.↵
 41.↵
 42.↵
 43.↵
 44.↵
 45.↵
 46.↵
 47.↵
 48.↵
 49.↵
 50.↵
 51.↵
 52.↵
 53.↵
 54.↵
 55.↵
 56.↵
 57.↵
 58.↵
 59.↵
 60.↵
 61.↵
 62.↵
 63.↵
 64.↵
 65.↵
 66.↵
 67.↵
 68.↵
 69.↵
 70.↵
 71.↵
 72.↵
 73.↵
 74.↵
 75.↵
 76.↵
 77.↵
 78.↵