Journal of Robotics, Networking and Artificial Life, Vol. 4, No. 3 (December 2017) 183–190___________________________________________________________________________________________________________

Force and Motion Analysis of larval zebrafish (Danio rerio)using a body dynamics model
Naohisa Mukaidani

Department of System Cybernetics, Graduate School of Engineering, Hiroshima University, 1-4-1 Kagamiyama,
Higashi-Hiroshima, Hiroshima, 739-8527, Japan
Zu SohDepartment of System Cybernetics, Institute of Engineering, Hiroshima University, 1-4-1 Kagamiyama, HigashiHiroshima, Hiroshima, 739-8527, Japan
Shinichi Higashijima

National Institutes of Natural Sciences, Okazaki Institute for Integrative Bioscience, National Institute for Physiological
Sciences, Okazaki, Aichi 444-8787, Japan
Toshio Tsuji

Department of System Cybernetics, Institute of Engineering, Hiroshima University, 1-4-1 Kagamiyama, HigashiHiroshima, Hiroshima, 739-8527, Japan
E-mail: mukaidani@bsys.hiroshima-u.ac.jp, sozu@bsys.hiroshima-u.ac.jp, shigashi@nips.ac.jp, tsuji@bsys.hiroshimau.ac.jp
www.bsys.hiroshima-u.ac.jp

Abstract
This study proposes a method of body dynamics analysis for larval zebrafish incorporating a viscoelastic bodymodel and a fluidic environment model to support the study of development mechanisms in motion generation. The
results showed that the estimated fluid drag coefficients enabled the body dynamics model to approximate the pathsof actual larvae with an error level of 0.76±0.74 [%] to the total body length.
Keywords: zebrafish, larval zebrafish, dynamics model, drag coefficientsystematically analyzed in the 1980s (Webb, 1984;
Webb and Weihs, 1986; Fuiman and Webb, 1988).
Subsequent development of computer vision techniques
enabled detailed quantification of larval motion, andfurther analysis helped to elucidate the ontogeny of finfunction (Danos and Lauder, 2008), the relationships
between drag force and fin/body morphologicaldevelopment from the larval stage to the adult stage(McHenry and Lauder, 2006) and the effects ofhydrodynamics on locomotor development (Danos,2012). Locomotor development was also analyzed in thecontext of larval body waves (Müller and Leeuwen,
2004), and the function of pectoral fins at the larval

1. Introduction
Zebrafish (Danio rerio), which measure about 3 – 40[mm] in length, have been widely used in the field of
biology to study genetic diseases and the brain structureof vertebrates. Their transparency at the embryonic andlarval stages of development facilitates in vivo
observation of their internal organs and neural circuits,
making them highly useful as a model organism for the
study of phenomena observed during the stages ofgrowth. For this reason, the behaviors of larval fish havebeen extensively investigated. For example, changes inswimming distance and velocity during ontogeny were

Copyright © 2017, the Authors. Published by Atlantis Press.
This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
183

Journal of Robotics, Networking and Artificial Life, Vol. 4, No. 3 (December 2017) 183–190___________________________________________________________________________________________________________

a background, this paper proposes a body dynamics
model incorporating the considerations of body form
and viscoelasticity to support the study ofdevelopmental mechanisms in the movement of larvalzebrafish. An algorithm for estimation of the externaldrag forces acting on such fish is also presented. The
proposed model enables simulation of target phenomena
in uniform virtual environments and with individual
conditions, and also allows analysis for the interactionof larval zebrafish morphology and swimming
movement.

stage was investigated by comparing the movements ofwild larval fish with those of others lacking pectoral fins
as a result of genetic manipulation (Green et al., 2011).
Based on the extensive body of behavioral analysis
data obtained from such studies, the mechanisms
supporting coordination between neural circuits and the
body for appropriate motor function have recently
attracted attention. For example, Budick and O’Malley
classified the behavior of larvae into three repertories(Budick, 2000), and Throsen et al. analyzed fin-axis
coordination (Throsen et al., 2004) to support
discussion of neural control. In addition, Higashijima et
al. succeed in indirect measurement of neural activities
in both the motoneurons and spinal interneurons during
escape behaviors using a calcium indicator calledcameleon (Higashijima, 2003).
The experimental studies described above succeededin elucidating a number of important mechanisms. The
clear-cut experimental results obtained have contributed
to enormous efforts regarding the control of complex
experimental conditions and static analysis to eliminate
the influences of individual differences in bodymorphology, learning ability and stages of development.
However, the effects of individual differences andexperimental conditions cannot be completely
eliminated in experimental approaches. In addition, it istechnically difficult to measure the forces produced bytiny larvae despite the importance of such informationin clarifying the neural control of muscles andswimming behavior. In light of such limitations,computer simulation can be seen as a potentiallyeffective approach. For example, Katsumata et al.
constructed a body dynamics model and carried out
fluid simulation for free-swimming larval zebrafish toevaluate the properties of fluid flow around them(Katsumata et al., 2009). This zebrafish larva modelallowed precise description of fluid behavior.
Considering the needs of biological study, however,both fluid dynamics and local body forces generated bylarvae are important in discussing the relationships
between ontogeny and nerve control development. In
this context, the development of a larva model that canbe used to calculate forces generated from the body andthose acting on the environment based on motion
recorded from actual fish can be seen as effective foranalyzing behavioral aspects of fish during ontogeny,
such as the relationships between drag force andmorphological development as analyzed by McHenryand Lauder (McHenry and Lauder, 2006). Against such

2. Body dynamics model
As larval zebrafish measure only about 3 [mm] inlength, the forces exerted between them and theirenvironments are difficult to measure. In this study, abody dynamics model was constructed based on thesizes and shapes of each part of these fish. A novel
method of estimating drag coefficients between larvalzebrafish and their environment as well as the driving
torque of each joint based on dynamic analysis
techniques was also introduced. Fig. 1 shows a zebrafish
larva at protruding-mouth stage (72 [h]) and a 3D
morphology of the body dynamics model configured
from this image. The proposed model is a rigid-link typecomposed of five parts: a head part consisting of one
link, left and right pectoral fin parts with 10 links each,a trunk part with 10 links, and a tail fin part with 10
links. The links are interconnected by joints allowing
rotation about the z-axis. The rotational motion of eachlink and the translational and rotational motion of a𝑇
representative point 𝑥g = �𝑥𝑔𝑔 , 𝜃𝑔 � of the model can becalculated from the following equation based on the
Newton-Euler method:
I(𝒒)𝒒̈ + h(𝒒̇ , 𝒒) + g(𝒒)= 𝝉 − 𝒌(𝒒 − 𝒒0 ) − c𝒒̇ + ∑𝑟,𝑑,𝑗 𝑱T𝑟𝑟𝑟 𝑭𝑟𝑟𝑟
I (𝒒)�𝒙 ̈ � + h (𝒒̇ , 𝒒) + g (𝒒) = ∑
𝑭
𝑔

𝑔

𝑔

𝑔

𝑟,𝑑,𝑗

𝑟𝑟𝑟

(1)
(2)

where I(q) and I𝑔 (𝒒 ) represent the inertia matrix,h(𝒒̇ , 𝒒) and h𝑔 (𝒒̇ , 𝒒) are the centrifugal force and the
Coriolis force, g(𝒒) and g 𝑔 (𝒒) represent the force ofgravity, 𝒒 = (𝒒𝒓𝒓 ) ∈ 𝐑30 is the rotational angle vectorof the joints, 𝒒0 is the equilibrium angle of the stiffness
component at a joint, 𝝉 = �𝝉𝒓𝒓 � ∈ 𝐑30 is a vector ofrotational driving force, k is the stiffness coefficient184

Journal of Robotics, Networking and Artificial Life, Vol. 4, No. 3 (December 2017) 183–190___________________________________________________________________________________________________________

Fig. 2. Viscoelasticity of the body dynamics model. The
white circles indicate inter-link joints, and qri denotes the angle
between the links. A pair of damper and spring units is
attached to each link in parallel to represent the viscoelastic
characteristics of muscles on the left and right sides. As only
two-dimensional motion was considered in the model for
simplification, the pairs of muscles on the dorsal and ventral
parts were combined.

centroid of a link as shown in the following equation(Mochizuki and Ichikawa, 2010):
Fig. 1. Shape of a larval zebrafish and the body dynamics
model. (a) shows a zebrafish larva at protruding-mouth stage(72 [h]) photographed by Kimmel et al. (Kimmel et al. 1995).
Thin transparent fins surround the dorsal and ventral parts ofthe trunk. (b) shows the configuration of the body dynamics
model approximating the image of the larva. The model’s
parameters were determined using this approximated shape.
Position x1 is the centroid of the head link, and x2 is the joint
connecting the head link and the trunk. The time traces ofthese two points were used to determine the drag coefficientsof a fluid environment.

𝑭𝑟𝑟𝑟 = 𝐶𝑟𝑟𝑟 𝜌𝑆𝑟𝑟𝑟 �𝒗𝑟𝑟𝑟 �𝒗𝑟𝑟𝑟

(3)

where 𝐶𝑟𝑟𝑟 is the drag coefficient, 𝜌 is the fluid density
and 𝑆𝑟𝑟𝑟 is the surface area.
As drag force significantly affects the model’s
motion, 𝐶𝑟𝑟𝑟 needs to be estimated. However, if
different drag coefficients were assigned to each of partr, direction d and link j, the number of coefficients to be
estimated would be as many as 82. Accordingly, it was
assumed for simplification that the drag coefficients forthe rostrocaudal direction of the pectoral fins and thetail fin (𝐶CA𝑗 , 𝐶RA𝑗 , 𝐶LA𝑗 ) were 0 because the fins are
negligibly thin compared to their total surface area. It
was also assumed that the drag coefficients for the same
part r were identical. As a result, the friction force
acting on each link can be determined via estimation of
6 drag coefficients.

(Fig. 2), c is the viscous damping coefficient (Fig. 2),
𝑱𝑟𝑟𝑟 is the Jacobian matrix, and 𝑭𝑟𝑟𝑟 is the friction force
exerted between the center of gravity of a link and theexternal environment. In addition, the subscript𝑟 ∈ (H, B, C, R, L) denotes parts of the model, with H
representing the head, B the trunk, C the tail fin, R theright pectoral fin and L the left pectoral fin. The
subscript 𝑑 ∈ (A, N) denotes direction, with A
representing the rostrocaudal direction of the fish and N
the normal direction in the rostrocaudal direction. The
subscript j denotes links sequentially numbered as
𝑗 = 1, 2, ⋯ , 10 from that connecting with the head part.
Although the parameters on the left of the motion
equation can be determined based on the measured
shapes and weights of actual larval zebrafish, drag force
is difficult to determine due to the size and complex
morphology of the fish. In this study, the drag force
𝑭𝑟𝑟𝑟 was approximated as a viscous resistance value
proportional to the square of the velocity 𝒗𝑟𝑟𝑟 at the

3. Drag coefficient estimation algorithm
As Eqs. (1) and (2) describe the relationship between
the motion of larval fish and related drag coefficients,
these coefficients can be estimated via dynamic analysis
by solving Eqs. (1) and (2). However, direct solution ofdrag coefficients from these equations for each timesample can yield large errors caused by video analysis
noise. Accordingly, the coefficients were determined inthis study by minimizing the following evaluation
function for robustness:

185

Journal of Robotics, Networking and Artificial Life, Vol. 4, No. 3 (December 2017) 183–190___________________________________________________________________________________________________________

Fig. 3. Four-step flow of the algorithm used to estimate thedrag coefficient of the environmental fluid. As dragcoefficients are related to the paths of the body dynamics
model, the proposed algorithm minimizes the difference
between the paths of the actual larva and the body dynamics
model. To avoid local minima, the algorithm sequentiallyperforms global search, local search and generalized reduced
gradient search where the drag coefficient sets with thesmallest errors were chosen at each step.
1

𝐸 = ∑𝑡 ∑𝑃𝑝=1�𝑥𝑝𝑅 (𝑡) − 𝑥𝑝 (𝑡)�𝑃

Fig. 4. Video image analysis procedure with recording
conducted at 1000 [fps]. (a) As the thinness of pectoral finsmakes it difficult to perform image analysis, the trunk and finsare manually separated and lines are drawn on the blurred
images of fins for interpolation. (b) Using the processed
images, joint bending angles and positions are calculated using
Wriggle Tracker software (Library Inc. Tokyo). (c) The
outcomes of motion calculation are input to the body
dynamics model, and the resulting paths are calculated basedon motion equations (Eqs. (1) and (2)) .

(4)

where P is the number of comparison points, 𝑥𝑝𝑅
represents positions on an actual fish, and 𝑥𝑝 represents

the corresponding positions on the body dynamics
model as shown in Fig. 1. This evaluation function
represents the error of paths between an actual fish andthe model as simulated using a set of drag coefficients.
The proposed algorithm minimizes this evaluation
function in M steps as shown in Fig. 3.
In Step 1, the movement of larval zebrafish isrecorded using a high-speed camera and video image
analysis is conducted (see Fig. 4). First, the tail andpectoral fins are separated as shown in Fig. 4 (a). Image
processing steps including binarization, noise reduction
and thinning are then performed for each frame of all
the separated images to extract the centroid line of thetrunk and the pectoral fins. This line is divided by points
corresponding to those of the model's joints, and thetime-dependent joint angles 𝑞𝑟𝑟 between the divided
points are calculated. As the measurement results
obtained from actual fish include noise, the joint angles
𝑞𝑟𝑟 are filtered using an n-degree low-pass Butterworth
filter. Finally, the joint angles are input to the model asshown in Fig. 4 (c).
In Step 2, a global search of the drag coefficients is
performed. First, a combination of drag coefficients is
𝑎1
𝑎1
𝑎1
𝑎1
𝑎1
𝑎1
selected from a set (𝐶HA1, 𝐶HN1, 𝐶BA1, 𝐶BN1, 𝐶RN1, 𝐶CN1)

defined using the equation below and substituted into
Eq. (3).
𝑎

𝑎

𝑎

𝑎

𝑎

𝑎

1
1
1
1
1
1
�𝐶HA1
, 𝐶HN1, 𝐶BA1, 𝐶BN1, 𝐶RN1, 𝐶CN1�1
1
1
1, 𝑁4 ∆𝐶BN1,∈ {𝑁1 ∆𝐶HA1 , 𝑁2 ∆𝐶HN1 , 𝑁3 ∆𝐶BA1
1
1
𝑁5 ∆𝐶RN1 , 𝑁6 ∆𝐶CN1 |𝑁1 , 𝑁2 , 𝑁3 , 𝑁4 , 𝑁5 , 𝑁61
= 0,1, ⋯ , 𝑁max}
(5)
1
where ∆𝐶𝑟𝑟1 is the step size for the global search
1yielding the 𝑎1 = 1,2, ⋯ , (𝑁max+ 1)6 variation of dragcoefficient sets, and the superscript 1 indicates the step
size used in the first optimization step. Data on the
motion of actual larvae, including the joint angle andjoint angular velocity of each link, are then substituted
into Eqs. (1) and (2), and the paths are calculated byrepeatedly determining the forward and inverse
dynamics for each sample time. This algorithm is
described below.
(i) The initial position and initial posture of arepresentative point on a larval zebrafish 𝒙𝑔 (0) arecalculated from image analysis results. Other initialvalues such as translational velocity 𝒙𝑔̇ (0), joint angle𝒒̇ (0), angular state velocity 𝒒0̇ (0) and angularacceleration 𝒒0̈ (0) are also calculated and substituted
into Eqs. (1) and (2).

186

Journal of Robotics, Networking and Artificial Life, Vol. 4, No. 3 (December 2017) 183–190___________________________________________________________________________________________________________(ii) The translational velocity 𝒗𝑟𝑟𝑟 (𝑡) of each link isderived from the translational velocity 𝒙𝑔̇ (𝑡) of arepresentative point and the angular velocity 𝒒̇ (𝑡) of
each joint. The translational velocity of each link issubstituted into Eq. (3) along with the selectedcombination of drag coefficients to determine thefriction force 𝑭𝑟𝑟𝑟 (𝑡).
(iii) The driving torque 𝝉(𝑡) and the acceleration 𝒙̈ (𝑡)are then derived from Eqs. (1) and (2) based on thesubstitution of 𝑭𝑟𝑟𝑟 (𝑡), 𝒒(𝑡), 𝒒̇ (𝑡) and 𝒒̈ (𝑡).
(iv) The position 𝒙(𝑡) of the representative point on
the model is derived via second-order integration for theacceleration 𝒙̈ (𝑡). By repeating steps (ii) through (iv)for each time step, the path is calculated for the selecteddrag coefficients. This path is then used to evaluate thecoefficients by comparing it to the path of the fish using
the evaluation function (Eq. (4)). After evaluation of all
combinations of drag coefficients given by Eq. (5), a set
𝑏1
𝑏1
𝑏1
𝑏1
𝑏1
of drag coefficients (𝐶HA1, 𝐶HN1, 𝐶BA1, 𝐶BN1, 𝐶RN1,𝑏1) with smaller 𝑑1 [%] values for path error E are
𝐶CN1
selected for a further search procedure in which1+ 1)6 𝑑1 /100.
𝑏1 = 1,2, ⋯ , (𝑁max
In Step 3, a local search is performed on theneighbors of the drag coefficient extracted in theprevious step. First, a set of combinations of drag
𝑎𝑚
𝑎𝑚
𝑎𝑚
𝑎𝑚
𝑎𝑚
𝑎𝑚
coefficients (𝐶HA1, 𝐶HN1, 𝐶BA1, 𝐶BN1, 𝐶RN1, 𝐶CN1) isselected. The search region for these coefficients isdefined by the following equation:
𝑎

𝑎

𝑎

𝑎

𝑎

min
In Step 4, a combination of drag coefficients (𝐶HA1
,
with
the
smallest
evaluation value E is searched using the generalized
reduced gradient method (Lasdon, 1974). The initialvalue for this method is selected from the drag𝑏2
𝑏2𝑏2
coefficients extracted in step 3, i.e., (𝐶HA1, 𝐶HN1, 𝐶BA1,𝑏2
𝑏2
𝑏2, 𝐶𝑅𝑁1, 𝐶𝐶𝐶1).
𝐶𝐵𝐵1With appropriate drag coefficients obtained using themethod described above, it is possible to determine thelocal force generated by the fish as the torque 𝝉 in Eq.
(1) driving each link.
min
min
min
min
min
𝐶HN1
, 𝐶BA1, 𝐶BN1, 𝐶RN1, 𝐶CN1)

4. Experiments
To verify the efficacy of the proposed method, dragcoefficient estimation was carried out using video
images of swimming larvae, and the effect of pectoral
fins on their swimming behavior was then explored.
4.1. Experimental condition
The experimental conditions configured for the dragcoefficient estimation algorithm (Fig. 3) are shown
below.
4.1.1. Video image analysis conditions
・Frame rate: 𝑓𝑠 = 1000 [fps]・Measurement duration: t m = 0.27 [s]・Low-pass filter: 𝑓𝑙𝑙𝑙 = 70 [Hz], n = 2
・The pectoral fin images used in the experiment weremanually interpolated for every frame due to an
insufficient level of clarity for analysis. In addition, thepectoral fins and the trunk were separated via manual
processing by painting unrelated parts of each frame inwhite as a background color (see Fig. 4 (b)).
・ Software used for video image analysis: Wriggle
Tracker (Library Inc., Tokyo)

𝑎

2
2
2
2
2
2
�𝐶HA1
, 𝐶HN1, 𝐶BA1, 𝐶BN1, 𝐶RN1, 𝐶CN1
�
𝑏1
𝑏1
𝑏1
2
2
, 𝐶𝐵𝐵1∈ {𝐶𝐻𝐻1 + 𝑁1 𝛥𝐶𝐻𝐻1 , 𝐶𝐻𝐻1 + 𝑁2 𝛥𝐶𝐻𝐻1
𝑏1
𝑏1
2
2
+𝑁3 𝛥𝐶𝐵𝐵1 , 𝐶𝐵𝐵1 + 𝑁4 𝛥𝐶𝐵𝐵1 , 𝐶𝑅𝑅1
𝑏1
2
2
, 𝐶𝐶𝐶1+ 𝑁6 𝛥𝐶𝐶𝐶1|𝑁1 , 𝑁2 , 𝑁3 ,+𝑁5 𝛥𝐶𝑅𝑅1𝑚
𝑚
𝑁4 , 𝑁5 , 𝑁6 = −𝑁max , ⋯ ,0, ⋯ , 𝑁max}
(6)2
where Δ𝐶𝑟𝑟1
is the step size of the drag coefficient2
1(Δ𝐶𝑟𝑟1 > Δ𝐶𝑟𝑟1) . The selected combination ofcoefficients is substituted into Eq. (3). Here, 𝑎2
represents the number of drag coefficients in the step 2.
Dynamic analysis is performed as described in Step 2,
and the simulated path is compared with that of an
actual larva using the evaluation function (Eq. (4)).
After evaluation of all combinations of the dragcoefficients shown in Eq. (5), another set of drag
𝑏2
𝑏2
𝑏2
𝑏2
𝑏2
𝑏2
coefficients (𝐶HA1, 𝐶HN1, 𝐶BA1, 𝐶BN1, 𝐶RN1, 𝐶CN1) withsmaller 𝑑2 [%] values of E is extracted for the end ofstep 3.

4.1.2. Configuration for drag coefficient estimation
1
・Step size: ΔCHA1= 1.3 × 10−8 ,1
1
Δ𝐶HN1 = 6.3 × 10−7 , Δ𝐶BA1= 1.3 × 10−7 ,1
1−4
Δ𝐶BN1 = 1.3 × 10 , Δ𝐶RN1 = 1.3 × 10−2 ,12= 1.3 × 10−4 , ΔCHA1= 1.3 × 10−9 ,
Δ𝐶CN1
2
2
Δ𝐶HN1 = 0.63, Δ𝐶BA1 = 1.3 × 10−8 ,2
2= 1.3 × 10−3 ,Δ𝐶BN1 = 1.3 × 10−5 , Δ𝐶RN1
2
Δ𝐶CN1 = 1.3 × 10−512・Number of searches: 𝑁max
= 100, 𝑁max
=4

187

Journal of Robotics, Networking and Artificial Life, Vol. 4, No. 3 (December 2017) 183–190___________________________________________________________________________________________________________

Fig. 7. Video images of the body dynamics model after
optimization. When optimal drag coefficients were chosen, themodel accurately traced the paths observed from an actual
larva.

Fig. 5. Video images of wild larval zebrafish motion. The
pictures show temporal changes in posture and position forevery 0.03 seconds from when the larva started swimming
until it stopped.

Fig. 8. Video images of the body dynamics model after
optimization. When optimal drag coefficients were chosen, themodel accurately traced the paths observed from an actual
larva

4.1.4. Configuration reading evaluation function
・Number of comparison points in evaluation function(Eq. 4): 𝑁 = 2
・Position of fish for evaluation (Eq. 4):
・Centroid of head (𝑖 = 1),
・Point of joint connecting head and trunk (𝑖 = 2) (seeFig. 1)

Fig. 6. Video images of the body dynamics model before
optimization. The pictures demonstrate a typical example withrandomly chosen drag coefficients in which drag force
precluded body thrust.

4.2. Drag coefficient estimation
Drag coefficients were estimated based on the algorithm
presented in Section 2.1. Fig. 5 shows video images of
the wild larval zebrafish. Fig. 6 shows the simulated
movement of the body dynamics model before
optimization, and Fig. 7 shows that after. Fig. 8 shows
the average path error between a wild fish and the body
dynamics model as given by Eq. (4), and indicates a
gradual decrease in each step. After drag coefficient
optimization using the generalized reduced gradientmethod, the path error approached 0 at about 0.76±0.74[%] compared to the body length. Fig. 9 shows the pathof larval movement, with (a) and (b) indicating a
position of 70 [%] normalized by body length with the
tip of the head as 0 [%] and the end of the tail as 100

4.1.3. Dynamic analysis conditions
・Equilibrium angle (Eq. 1): 𝑞0 = 0・Density of water (Eq. 3): 𝜌 = 1.00 [g/cm3 ](at a water temperature of 20 [℃])・Larvalmodel
length:𝑙𝑏 = 2.78 [mm]・Mass density of larval model: 3.33 [g/cm3 ](Calculated based on reference of Avella et al. (2012))
・Form of link sets: as shown in Figs. 1 and 4
・Joint stiffness: 15.5 [kNm ⋅ rad/m2 ] (Long, 1998)
・Joint viscosity: 0.84 [kNms ⋅ rad/m2 ] (Long, 1998)

188

Journal of Robotics, Networking and Artificial Life, Vol. 4, No. 3 (December 2017) 183–190___________________________________________________________________________________________________________

Fig. 10. Resultant force of the trunk and fins. The broken line
shows the resultant force of the pectoral fin part, and the solidline shows that of the trunk and tail part. Overall, the force
generated from the pectoral fins is negligible compared to thatgenerated from the trunk, which is consistent with biologicalexperiment results [Green et al. 2011]. However, in period (a),the larvae initiated swimming and the model predicted that the
pectoral fins would generate force to push the fish backward.

Fig. 9. Comparison of paths obtained before and after
optimization. (a) shows a typical path of the body dynamics
model when drag coefficients were chosen randomly, andcorresponds to the images shown in Fig. 6. (b) demonstrates
the path of the body dynamics model given optimal dragcoefficients, and corresponds to the images shown in Fig. 7. It
is particularly notable that the model successfully traced the
backward swimming observed when the actual fish began
swimming.

The resultant force generated by the trunk and tail fin
resrespart (𝑭res
BC = 𝑭B +𝑭C ) and that of the pectoral fin part
res
res
res
(𝑭RL = 𝑭R +𝑭L ) were calculated. In Fig. 10, period(a) shows the time immediately after the larva started
swimming (0 – 0.045 [s]) during which it moved
backward, and period (b) shows the duration of forward
swimming (0.045 – 0.25 [s]). For period (b), theanalysis results showed that the average resultant force� resof the trunk and tail fin part was about 𝑭
RL =78 [𝜇N]
� resand that of the pectoral fin part was about 𝑭
BC =4068[ 𝜇N ]. As the resultant force of the pectoral fin part
accounts for only 2 [%] of the total force generated by
the body dynamics model, its influence on swimming
can be considered negligible. This result is consistent
with those of the experimental studies. Green et al.
assumed that pectoral fins affected the movement oflarval zebrafish and conducted motion analysis on
actual fish (Green et al., 2011). The results showed aminimal influence from pectoral fins on swimming in
the larval stage; rather, the fins are used to supply
oxygen to the gills at this stage of development. Thisoutcome verified the efficacy of the proposed model anddrag coefficient estimation algorithm, indicating their[%]. The solid line shows the path of the larvalzebrafish, and the dashed line shows the path of themodel before and after drag coefficient optimization.
Comparison of the optimized model and the larvalzebrafish shows that the former captured thecharacteristics of fish motion in which both the modeland the larval zebrafish started by swimming backward
then forward. These results indicate that the body
dynamics model with drag coefficients estimated using
the proposed algorithm successfully reproduces the
movement of wild larval zebrafish.
5.

Results and Discussion

Although the pectoral fins of adult zebrafish are known
to contribute to thrust (Webb, 1973), their function indevelopmental-stage larval zebrafish has not been fully
elucidated. In this study, the influence of pectoral finson swimming was examined by comparing the resultant
force generated from the pectoral fin part with those ofother parts. These resultant forces were calculated bysimply adding all force vectors acting on thecorresponding links via the following equation:𝑭res𝑟 = ∑𝑑𝑑 𝑭𝑟𝑟𝑟

(7)

189

Journal of Robotics, Networking and Artificial Life, Vol. 4, No. 3 (December 2017) 183–190___________________________________________________________________________________________________________

3. N. Danos and G. V. Lauder, The ontogeny of fin function
during routine turns in zebrafish Danio rerio, J. Exp. Biol.
210, (2008), 3374-3386.
4. L. Fuiman and P. W. Webb, Ontogeny of routine
swimming activity and performance in zebra danios(Teleostei: Cyprinidae), Anim. Behav. (1988), 36, 250261.
5. M. H. Green, R. K. Ho and M. E. Hale, Movement andfunction of the pectoral fins of the larval zebrafish (Daniorerio) during slow swimming, Exp. Bio. 214, (2011),
3111-3123.
6. S. Higashijima, Y. Hotta and H. Okamoto, Visualizationof cranial motor neurons in live transgenic zebrafish
expressing green fluorescent protein under the control of
the islet-1 promoter/enhancer, J. Neurosci, 20, (2000),
206–18.
7. Y. Katsumata, U. K. Muller and H. Liu, (2009).
Computation of Self-Propelled Swimming in Larva
Fishes, Bio. Sci. and Eng. 4, 1, 54-66.
8. C. B. Kimmel, W. W. Ballard, S. R. Kimmel, B. Ullmann
and T. F. Schilling, Stages of Embryonic Development of
the Zebrafish, Dev. Dynam. 203, (1995), 253-310.
9. L. S. Lasdon, R. L. Fox and M. W. Ranter, Nonlinear
optimization using the generalized reduced gradientmethod, Revue Francaise d’Automatique. Informatique etRecherche Operationelle. 3, (1974), 73-104.
10. J. H. Long Jr. Muscles, Elastic Energy, and the Dynamics
of Body Stiffness in Swimming Eels, Amer. Zool. 38,(1998), 771-792.
11. M. McHenry and G. Lauder, Ontogeny of form andfunction: locomotor morphology and drag in zebrafish(Danio rerio), J. Morphol. 267, (2006), 1099-1109.
12. O. Mochizuki and S. Ichikawa, Hydrodynamics forlearning from creatures (in Japanese), Yokendo, (2010).
13. U. K. Müller and J. L. van Leeuwen, Swimming of larvalzebrafish: ontogeny of body waves and implications forlocomotory development, J. Exp. Biol. 203, (2004), 853868.
14. D. H. Thorsen, J. J. Cassidy and M. E. Hale, Swimmingof larval zebrafish: fin-axis coordination and implication
for function and neural control, J. Exp. Biol. 207, (2004),
4175-4183.
15. P. W. Webb, Effects of partial caudal-fin amputation on
the kinematics and metabolic rate of underyearlingsockeye salmon (Oncorhynchus nerka) at steadyswimming speeds, J. exp. Biol. 59, (1973), 565-581.
16. P. W. Webb, Form and function in fish swimming.
Scient. Am., 251, (1984), 72-82.
17. P. W. Webb and D. Weihs, Functional locomotor
morphology of early life history stages of fishes, Trans.
Am. Fish. Soc. 115, (1986), 115-127.

suitability for the prediction of unknown dynamic
characteristics. In this context, focus was subsequently
placed on period (a), for which the average resultant
� resforce of the trunk and tail fin part was about 𝑭
BC =39[ 𝜇N ] and that of the pectoral fin part was about
� res
𝑭
RL =121 [𝜇N]. These results indicate that the trunk
and tail fin part have a smaller influence in period (a),allowing pectoral fin movement to contribute more toswimming. Indeed, the fish was observed to maneuver
its pectoral fins toward its head, thereby creating force
for backward swimming. Such pectoral fin movement
may thus be applied to generate force in order todetermine the direction of swimming. In this way, theproposed model can provide guidelines for biologicalexperiments.
6. Conclusions and future work
This paper proposed a body dynamics model for larvalzebrafish and reported on dynamic analysis of thespecies. The measured paths of wild larval fish werecompared with those obtained from swimming
simulation using the proposed model. The results
indicated that the model can be used to approximate
paths with a high accuracy level of 0.76±0.74 [%] in
relation to body length. It was also found to be capableof reproducing paths with characteristics similar to
those of wild larval fish as well as supporting estimationof how pectoral fins affect swimming.
In future work, the author plans to elucidate theswimming mechanism of larval fish with focus on the
motion generation mechanisms hidden behind the neuralcircuit. To this end, muscle and neural circuit models
will be constructed, and the mechanisms behind thedevelopment of swimming behavior in larval zebrafish
will be analyzed
Reference
1. M. A. Avella, A. Place, S.-J. Du, E. Williams, S. Silvi, Y.
Zohar and O. Carnevali, (2012). Lactobacillus rhamnosus
accelerates zebrafish backbone calcification and gonadal
differentiation through effects on the GnRH and IGF
systems, PLoS One 7, (2012), e45572
2. S. Budick and D. O’Malley, Locomotor repertoire of thelarval zebrafish: swimming, turning, and prey capture, J.
Exp. Biol. (2000), 203, 2565-2579.

190