A Procedural Method for Modeling the Purkinje Fibers of the Heart
TAKASHI IJIRI1, TAKASHI ASHIHARA2, TAKESHI YAMAGUCHI3, KENSHI TAKAYAMA1, TAKEO IGARASHI1,TATSUO SHIMADA3, TSUNETOYO NAMBA4, RYO HARAGUCHI5, KAZUO NAKAZAWA5
1

Department of Computer Science, the University of Tokyo
Department of Cardiovascular Medicine, Heart Rhythm Center, Shiga University of Medical Science3
Department of Health Science, School of Medicine, Oita University
4
Department of Medical Technology, Kagawa Prefectural College of Health Sciences
5
National Cardiovascular Center Research Institute2

away from existing branches as much as possible tocreate a uniform distribution, and 2) when branches
collide, we connect the colliding branches to construct a closed mesh structure. We designed a generating rule based on observations of the photograph of Purkinje fibers and manually specified
three terminal positions on a three-dimensional (3D)
heart model: those of the right bundle branch, theanterior fascicle, and the left posterior fascicle of theleft branch. Then, we grew fibers starting from eachof the three positions based on the specified generating rule.
We achieved to generate 3D models of Purkinje
fibers of which physical appearances closely resembled the real photograph. The generation takes
a few seconds. Variations of the Purkinje fibers
could be constructed easily by modifying the generating rules and parameters.

Abstract: The Purkinje fibers are located in the ventricular walls of the heart, just beneath the endocardium and conduct excitation from the right and left
bundle branches to the ventricular myocardium. Recently, anatomists succeeded in photographing the
Purkinje fibers of a sheep, which clearly showed themesh structure of the Purkinje fibers. In this study,
we present a technique for modeling the meshstructure of Purkinje fibers semiautomatically using
an extended L-system.
The L-system is a formal grammar that defines
the growth of a fractal structure by generating rules(or rewriting rules) and an initial structure. It was
originally formulated to describe the growth of plant
cells, and has subsequently been applied for various
purposes in computer graphics such as modelingplants, buildings, streets, and ornaments. For our
purpose, we extended the growth process of the Lsystem as follows: 1) each growing branch keeps
Key words: Purkinje fibers, L-System, heart simulation.
---INTRODUCTION
A three-dimensional (3D) virtual heart model is often
used for computer simulations and visualizations. Computer simulation is one way to understand the electrophysiological properties of the heart or to figure out the mechanisms of fatal arrhythmias [1-4]. Effective visualization of a 3D heart model is also a useful tool for education and communication between doctors and patients [1,5]. However, the creation of a 3D heart model is difficult
and time-consuming because the heart has intricate structures containing various tissues, such as the atrioventricular node, bundle of His, Purkinje fibers, and contractive
myocardium. Our goal was to facilitate this process byproviding effective modeling tools. In this study, we focused on the construction of Purkinje fibers.
The Purkinje fibers are part of the ventricular conduction system and were originally discovered by Tawara [6].
These tissues conduct excitation (electrical activation)
rapidly from the bundle of His to the ventricular myocardial tissue. The Purkinje fibers are located in the ventricular walls of the heart, just beneath the endocardium. Figure 1 is a PAS-stained stereomicrograph of a sheep heart
provided by Shimada et al. [7], which shows the closed
mesh structure of the Purkinje fibers. Since the Purkinje
fibers play a key role in the conduction system, many
simulation studies have taken them into account. In such
simulation-based systems, the fibers are usually constructed manually. In [2, 8, 9], authors drew a simplified

1mm
Fig. 1. A photograph of the Purkinje fibers of asheep [7].
two-dimensional (2D) shape of the Purkinje fibers manually on the flatten endocardium and then mapped the
2D shape onto a 3D heart model. However, it is still difficult to manually create detailed structures as shown in Fig.
1. It also takes time to generate various Purkinje fiber
patterns which are sometimes required to perform simulations using various models.
In this report, we present a semiautomatic system formodeling the mesh structure of the Purkinje fibers by applying the L-system. The user manually specifies the endocardial regions that the Purkinje fibers will cover, sev-

eral positions from which fibers will grow (e.g., the terminal positions of the right bundle branch, the left anterior fascicle, and the left posterior fascicle on the endocardium), and the parameters of the L-system. Then, thegrowth simulation of L-system generates 3D mesh structures that cover the specified regions. This process takes
only a few seconds and the user can generate various
Purkinje fiber models by simply changing the generation
parameters.

Branch segment

(a)

(b)

Fig. 3. Branch segments of the Purkinje fibers(a) and volume representation of the Purkinje
layer (b).

click

(d)

(a)

Fig. 2. A simple example of the L-system (compound leaf) [12]. Beginning with the initial structure(leftmost thin red line), the system generates acomplicated structure by applying two generating
rules (top-left) sequentially.
L-system: The L-system is a formal grammar originallyintroduced by Lindenmayer to formalize the developmentof multicellular plants [10] and was subsequently expanded to represent higher plants and complex branching
structures [11-13]. The framework of the L-system consists of an initial structure and rewriting rules (or generating rules). The essence of development is parallel replacement using the rewriting rules. Starting from theinitial structure, the L-system replaces each part of thecurrent structure by applying the rule sequentially. Figure
2 shows a simple example, the development of a compound leaf [12]. This includes two module types: theapices (thin red lines) and the internodes (thick blacklines). This example has two rewriting rules (Fig. 2, topleft): one replaces an apex with an internode, two lateralapices, and a top apex, while the other replaces an internode with a longer one. The initial structure is a single
apex. Using these simple rules, the system develops an
intricate branching structure over several replacing steps.
An interesting aspect of the system is that each replacement process corresponds to the growth of part of theplant. Therefore, the L-system is not only a heuristic
technique that creates fractal-like shapes, it is also a simulation of real-world plant growth. Karch et al. presented asimilar method to generate vessel tree patterns [14].
However, the L-system and the method of Karch et al.
[14] are designed for open tree structures and cannot beapplied directly for the closed mesh structures seen in
Purkinje fibers.
METHODS
In contrast to existing approaches that create Purkinje
fibers on flattened 2D endocardium [2, 8, 9], we generatethem directly in 3D space to achieve the detailed 3D
overlapping structures of Purkinje networks. Our systemloads a 3D polygon model of the heart and generates Purkinje fiber models just beneath the endocardium. We

(b)

(c)

(e)

Fig. 4. Overview of the modeling process. The
user specifies endocardial regions (b) on a 3D heartmodel (a) so as to create the Purkinje layer (c). The
user also specifies the initial points by clicking (d,e).

apex
segment(a)(b)

Fig. 5. A generating rule (a) and its handle (b).
The user can modify the growth directions of twochild apices by manipulating the control points.
represent the mesh structure of the Purkinje fibers with anumber of branch segments connected to each other (Fig.
3a). Each branch is represented as a curved polyline andconstructs cylindrical objects along the polyline for effective visualization. We also create a volume representation
to compute a distance field to guide the growth of the
Purkinje fibers. We prepare 128 × 128 × 128 volume datathat cover the entire heart model, and activate voxels justbeneath the endocardium where the Purkinje network is
located (we call this the Purkinje layer; Fig. 3b).
User interface: After loading a 3D heart model, theuser first specifies endocardium surface regions on theheart model. We provide a sketch-based interface thatallows the user to specify polygons by painting. The specified surfaces are highlighted in red (Fig. 4b). The system then generates the Purkinje layer by activating voxelsthat are inside the 3D heart model and close to the specified endocardium surface (we used 3r as a threshold,where r is the voxel width; Fig. 4c). Next, the user specifies initial points from which the Purkinje fibers willgrow by clicking on the heart model directly. The system
places four initial apexes radially at each specified position (Fig. 4d, 4e). Finally, the user sets the growth parameters of the generating rule by manipulating the control
points directly [15] (Fig. 5). When all of the parametersare set, the user pushes the “growth button” to start thegrowth simulation of the L-system.
The actual photograph of Purkinje fibers (Fig. 1) showsthat most of the branch segments are divided into twochild branches. Based on this observation, we provide a
simple generating rule in Fig. 5. This generating rule contains two module types, apices (blue arrow) and internodes (orange bar), and replaces one apex with an internode and two new child apices. The system allows theuser to modify directions of the two child apices by manipulating a handle (Fig. 5b). The user can also set branch
length lbra. We introduce a randomness to the length byusing Gaussian random number: lbra’ = lbra + random( μ,σ2). We set mean value μ = 0, and standard deviation σ2 =0.4 × lbra for generating the results in this paper.
Creation of Purkinje Network: When the growth button is pressed, our system runs a growth simulation tocreate a Purkinje network. The growth simulation uses aniterative algorithm similarly to the standard L-system.
Each iteration process consists of two steps. In the firststep, we randomly shuffle the order of all apices and insert them into a queue. (Note that there are 4 apices foreach user-specified initial point at the beginning.) In thesecond step, we grow all apices in the queue one by one.
Differently from the standard L-system in which allbranch segments have straight shapes (Fig. 2), we growan apex to be a curved branch segment so as to obtain auniform distribution (see the next paragraph). If a growing apex steps out of Purkinje layer or collides to an existing branch, we stop the growth of the apex. If not, we
add two new apices based on the generating rule (Fig. 5).
These new apices will grow in the next iteration. Theirdirections and lengths are determined by the userspecified parameters. We stop the growth iteration whenno new apices are generated.
Extended L-system: The branch segments of Purkinje
fibers are curved and distributed uniformly, as shown in
Fig. 1. To generate such mesh structures, we introduce
two extensions of L-system: one is for generating uniform
distributions and the other is for constructing closed meshstructures.
To generate uniform branch distributions, we grow andbend a new branch to keep away from all existing
branches. Our system approximates the curved shape of abranch with a simple polyline consisting of five line segments. During the growth process, the system first constructs a distance field from all existing branches and calculates the gradient of the distance field (Fig. 6a). We
compute only this distance field and the gradient field at
the voxel grids for fast calculations. The system thengrows a new branch that curves along the gradient. We
define the direction of each line segment of a branch as;(1)d=d original + w1 * d gradient|| d original + w1 * d gradient ||.

(1)

where doriginal is the direction of the previous segment,
dgradient is the gradient direction at the terminal point of
the previous segment, and w1 is a user-specified weight.
Figure 6b shows a simple example: a new branch is growing from bottom to top in a gradient field that slants fromleft to right. The branch curves gradually in the directionof the gradient. We also project the direction d parallel tothe endocardium surface so that the growing branch does

dgradient(a)doriginal(b)(c)

Fig. 6. The distance field of one branch (a) andthe growth of a branch in a gradient field (b, c).

Fig. 7. Excitation sequence on the ventricular
endocardium. The numbers are the times in milliseconds when stimuli occur at points [16].
not pop out from the thin Purkinje layer:

d ' = d − (d ⋅ n) n .

(2)
where n is the normal vector of the nearest polygon on
the endocardium surface. When a growing branch steps
out of the user-specified endocardial region, the systemstops its growth.
To generate a closed mesh structure, we add a simple
rule to the L-system. If a growing branch collides with anexisting branch, the system connects the growing branch
to the collided branch and stops the growth of the branch.
Since thin branches rarely collide with each other in 3D
space, we use a threshold k. If the distance between growing and existing branches is less than k, we detect a collision. We set k = 0.2 × lbra.
RESULTS AND DISCUSSION
Based on a photograph of actual Purkinje fibers, we designed a generating rule (Fig. 5a) and set branch lengthlbra = 1.8 mm. The branch length is ~10-times the actual
one since the Purkinje network has self-similarity as
shown in Fig. 1. The conduction velocity in the Purkinje
fibers is ~2 m/s and thus the activation time along thebranch is estimated as ~0.9 ms, which is markedly shorter
than the action potential rising phase (several milliseconds). Therefore, we strongly believe that this degradation of the Purkinje network scale does not alter the physiological role in the simulated 3D ventricles. Next, we
specified 3 initial growth points based on a map of excitation sequence on the endocardium that was modeled
based on measured data (Fig. 7) [16]. These 3 initial
growth points (Fig. 8a) correspond to the 3 earliest stimuli regions (0 ms regions in Fig. 7); distal ends of the rightbundle branch, the left anterior fascicle, and the left posterior fascicle. Figure 8b shows a generated Purkinje model, demonstrating that our system can capture the detailed
mesh structure. The growth simulation takes less than 10

(a)

(b)

Fig. 8. The 3D heart model (a) and 3D Purkinje fibers created by our system (b). In (a), the endocardial
regions are highlighted in red and the initial points of growth are in yellow. The resulting model, whose physical appearance closely resembles that of the actual Purkinje fibers, is generated by our system.
% of branches

20

real photo10

our model

5
0
0~2
2~4
4~6
6~8
8~10
10~12
12~14
14~16
16~18
18~20
20~22
22~24
24~26
26~28
28~30

1mm

1mm

15

Length of branch (10μm)

(a)

(b)

(c)

(d)

Fig. 9. An evaluation of the resulting model. A chart (d) compares a branch length distribution of real Purkinje fibers in the photograph and that of our resulting model.
In contrast to existing methods that require the user toseconds.
model the detailed mesh structures manually, our system
To evaluate the quality of Purkinje fiber model createdallows the user to design variants of 3D Purkinje fiber
by our method, we compared a distribution of branchmodels using simple interactions. Since the growth simulengths of our model and that of the real photograph inlation using the extended L-system is rapid, the user can
Fig. 1. First, we manually binarized the photograph (Fig.
easily examine and tune parameters by trial and error. The9a) and extracted mesh structures by skeletonizing theresulting Purkinje fiber models have detailed mesh strucbinary image (Fig. 9b). Next, we generated real-scaletures and their physical appearances closely resemble the
Purkinje fibers on a plane surface by our method (Fig. 9c).
photograph of the Purkinje fibers.
We set lbra = 0.18 mm. We then compared the distributionof branch lengths (Fig. 9d), and found that our created
CONCLUSIONS
model has a similar length distribution to the real Purkinje fibers. In both cases (Figs. 9b and 9c), the mean
We present a semiautomatic system for creating 3D Purlength of branches is ~0.13 mm. In the generation process
kinje fiber models and introduce two extensions to the Lof Fig. 9c, we observed that ~45% apices collided to exsystem: collision avoidance and closed mesh structure
isting branches and stopped their growths. This causes
creation. The user specifies the endocardial regions, inithat the number of shorter branches are more than that oftial points of growth, and generating rules of the L-system.
longer branches. Figure 10 shows the growth iteration
Then, our extended L-system performs the growth simuprocesses of the extended L-system on a plane surface.
lation that constructs the 3D mesh structure of Purkinje
We can also generate variations of Purkinje fiber models
fibers in a few seconds.
by modifying the parameters of the L-system, as shown in
One limitation is that our system requires users to tune
Fig. 11.
the parameters based on their observations and experiWe can also perform a simple simulation of excitation
ments. We plan to extend this system so that rules andconduction on Purkinje fibers. First, we activate branches
parameters can be estimated from anatomical photoconnected to the initial points. Then, we iteratively actigraphs. Another future project is to apply our method tovate branches connected to the activated branches in eachother targets such as airways and blood vessels.
simulation step. Figure 12 illustrates the excitation conduction along the Purkinje fibers. We observe an excitaThis work was supported in part by Adobe Systems Inc.
tion conduction pattern similar to the measured modeland JSPS Research Fellowship.
(Fig. 7).

Growth direction

Initial point
Fig. 10. Growth process of the extended L-system on a plane surface. These panels show the 1st, 3rd,
6th, 9th, and 12th iterations of the growth process from left to right. Starting from the center point, the Purkinje
fibers grow to fill the pink region.

GrowthDirection

w1 = 0.75
lseg= 23μm

w1 = 0lseg= 23μm

(a)

(b)

w1 = 0.75
lseg= 23μm(c)

w1 = 0.75
lseg= 23μm

(d)Fig. 11. Variations of Purkinje fiber models and their growth parameters. In example (a), we turn off thecollision avoidance framework by setting w1 = 0.

LV

RV

Fig. 12. Excitation conduction along the ventricular Purkinje fibers. The excitation conduction started
from the initial points of growth, and we activated the connected branches iteratively. These panels show the
2nd, 6th, 10th, 14th, 18th, and 22nd iterations.
REFERENCES
1. Nakazawa K, Suzuki T, Ashihara T, Inagaki M, Namba T, Ikeda T,
Suzuki R: Computational analysis and visualization of spiral wave
reentry in a virtual heart model. In: Clinical Application of Computational Mechanics to the Cardiovascular System, ed. Yamaguchi
T, Springer, New York, 2000
2. Berenfeld O and Jalife J: Purkinje-muscle reentry as a mechanism of polymorphic ventricular arrhythmias in a 3-dimensionalmodel of the ventricles. Circ Res 82: 1063–1077, 1998.
3. Smith N and Hunter P: Giving form to the function of the heart:
embedding cellular models in an anatomical framework. JJP 54:
541–544, 2004
4. Usyk TP, LeGrice IJ, and McCulloch AD: Computational model ofthree-dimensional cardiac electromechanics. Comput Visual Sci
4: 249–257, 2002
5. Igarashi T, Ashihara T, Yao T, Nagata S, Takada M, Sakachi H,
Suzuki T, Nakazawa K: Examination of various pen-based input
methods for electronic medical recording systems. In Proc. of thest
21 Joint Conference on Medical Informatics, 2001
6. Tawara S: Das Reirzleirungssysrem des Saugertierherzens.
Fischer, JENA, 1906
7. Shimada T, Ushiki T, and Fujita T: Purkinje fibers of the heart.
Shinyaku to Chiryou 42: 11–13, 1992 (in Japanese)
8. Clements CJ and Vigmond EJ: Construction of a cardiac conduction system subject to extracellular stimulation. In Proc. of theth2005 IEEE Engineering in Medicine and Biology 27 Annual Conference Shanghai, 4235–4238, 2005
9. Pollard AE and Barr RC: The construction of an anatomically
based model of the human ventricular conduction system. IEEE
Trans Biomed Eng 37: 1173–1185, 1990
10. Lindenmayer A: Mathematical models for cellular interactions indevelopment, I & II. J. Theor Biol 18: 280–315, 1968
11. Prusinkiewicz P and Lindenmayer A: The Algorithmic Beauty ofPlants. Springer-Verlag, New York, 1990
12. Prusinkiewicz P, Hammel M, Hanan J, and Měch R: L-systems:ndfrom the theory to visual models of plants. In Proc. of the 2
CSIRO Symposium on Computational Challenges in Life
Sciences, 1996
13. Měch R, and Prusinkiewicz P: Visual models of plants interactingwith their environment. ACM SIGGRAPH 96, 397–410,1996
14. Karch R, Neumann F, Neumann M, and Schreiner W: A threedimensional model for arterial tree representation, generated byconstrained constructive optimization. Comput Biol Med 29: 19–
38, 1999
15. Ijiri T, Owada S, and Igarashi T: The Sketch L-System: global
control of tree modeling using free-form strokes. In Proc. of Smart
Graphics 2006, 138–146, 2006
16. Nakazawa K, Ashihara T, Suzuki T, Yao T, and Namba T: Cardiac
model based on active membrane model: spiral wave reentry asa mechanism of lethal arrhythmia. In: Physiome for the Cariovascular System, ed. Okamoto Y. Morikita, Tokyo, 2003(in Japanese)