We developed a 2D heat flux model to elucidate routes and rates of heat transfer within bigeye tuna Thunnus obesus Lowe 1839 in both steady-state and time-dependent settings. In modeling the former situation, we adjusted the efficiencies of heat conservation in the red and the white muscle so as to make the output of the model agree as closely as possible with observed cross-sectional isotherms. In modeling the latter situation, we applied the heat exchanger efficiencies from the steady-state model to predict the distribution of temperature and heat fluxes in bigeye tuna during their extensive daily vertical excursions. The simulations yielded a close match to the data recorded in free-swimming fish and strongly point to the importance of the heat-producing and heat-conserving properties of the white muscle. The best correspondence between model output and observed data was obtained when the countercurrent heat exchangers in the blood flow pathways to the red and white muscle retained 99% and 96% (respectively) of the heat produced in these tissues. Our model confirms that the ability of bigeye tuna to maintain elevated muscle temperatures during their extensive daily vertical movements depends on their ability to rapidly modulate heating and cooling rates. This study shows that the differential cooling and heating rates could be fully accounted for by a mechanism where blood flow to the swimming muscles is either exclusively through the heat exchangers or completely shunted around them, depending on the ambient temperature relative to the body temperature. Our results therefore strongly suggest that such a mechanism is involved in the extensive physiological thermoregulatory abilities of endothermic bigeye tuna.

Tuna are highly specialized predatory fish with a unique assembly of anatomical, physiological and biochemical adaptations. Numerous studies have investigated these specializations such as fusiform body shape, thunniform swimming mode, capacity to conserve metabolic heat in red muscle, elevated metabolic and physiological rate functions, and a frequency-modulated cardiac output (Carey et al., 1984; Stevens and McLeese, 1984; Graham and Dickson, 2000; Lowe et al., 2000; Altringham and Shadwick, 2001; Brill and Bushnell, 2001; Gunn and Block, 2001; Katz et al., 2001; Sepulveda et al., 2003; Blank et al., 2004). In particular, the ability to conserve metabolically produced heat due to centralized red muscle and the presence of well-developed vascular countercurrent heat exchangers (also known as `retes') has been thoroughly investigated (e.g. Carey et al.,1971; Mathieu-Costello et al.,1992; Block et al.,1993; Holland and Sibert,1993; Dewar et al.,1994; Graham and Dickson,2001; Korsmeyer and Dewar,2001). Most attention has been focused on the heat exchangers associated with the red muscle fiber portions of the myotomes (herein referred to as `red muscle'). However, the white muscle fiber portions of the myotomes(herein referred to as `white muscle') account for approximately 45–55%of total body mass (Graham et al.,1983; Freund,1999) and the potential significance of the capacity of white muscle to produce and conserve heat remains largely unexplored.

A growing body of evidence confirms the ability of tuna to modulate muscle temperature by physiological thermoregulation or by a combination of physiological adjustments and behavior. This allows these fishes to appropriately adjust rates of body temperature change in the face of rapid ambient temperature changes induced by their extensive daily vertical excursions (Dewar et al.,1994; Graham and Dickson,2001; Malte et al.,2007). It has been demonstrated that free-ranging bigeye tuna can change their whole-body thermal conductivity (as expressed by heating and cooling rates) by a factor of 100–1000 during vertical excursions(Holland and Sibert, 1993; Musyl et al., 2003). Modeling these observations in a two-compartment lumped parameter model further defined the ability of bigeye tuna to maintain elevated red muscle temperatures through behavioral and physiological means(Malte et al., 2007). The simple non-dimensional model described by Malte and colleagues(Malte et al., 2007) did not,however, distinguish between the different tissues in the tuna (the red muscle temperature referred to is just an overall elevated body temperature). Hence the model cannot reveal where in the body the heat is produced or conserved.

Carey and Teal (Carey and Teal,1966) and Carey (Carey,1973) confirmed the presence of an advanced countercurrent heat exchanger in red muscle when examining bigeye tuna, but they also describe,`...bands of alternate arteries and veins running into the light [i.e. white]muscle. These bands are widest and most numerous near the lateral vessels,where they may be 20 vessels wide, but decrease to triads of artery–vein–artery toward the ventral and dorsal midlines'. The presence of a simple countercurrent heat exchanger in the white muscle reported in the early literature clearly provides evidence for the potential role of the white muscle in heat conservation and thermoregulation. We therefore decided to investigate the heat flux in different regions in the bigeye tuna. Because it is cumbersome to obtain direct measurements of the body temperature of tuna in motion (Dewar et al., 1994; Holland and Sibert, 1993), and because no experiment has so far elucidated the transectional heat fluxes of any tuna, we decided to address these questions through mathematical modeling of heat flux. The questions we specifically addressed were: (1) does the white muscle contribute to the conservation of heat? (2) If so, is this conservation even noticeable? By developing a 2D model, we were able to simulate the heat flux and accurately predict temperature measurements made in free-swimming tuna, thereby acquiring a more profound understanding of the precise mechanisms responsible for the thermoregulatory abilities in this species. Our approach also enabled us to investigate a variety of features, including qualifying and quantifying the significance of the white muscle with respect to heat transfer and thermoregulation. We also determined the routes of heat loss, the possible extent of the modulation of heat transfer rates, and the role of body size.

To explore the behavioral and physiological thermoregulatory capacities of bigeye tuna, two different models and related sets of data were required. (1)A steady-state model simulating the transectional heat flux through the tuna's heat-producing tissues and heat-conserving modes. (2) A non-steady-state time-dependent model simulating the events associated with bigeye tuna's daily vertical movements. For the former we used data on isotherms in a cross-section of a bigeye tuna from Carey and Teal(Carey and Teal, 1966). For the latter we employed data on body temperature, water temperature and depth in free-ranging bigeye tuna from Musyl et al.(Musyl et al., 2003). In brief, Carey and Teal (Carey and Teal,1966) caught bigeye tuna on long-line fishing gear during the winter of 1965–66 in the western North Atlantic and the eastern South Pacific. Once on board, the fish were sedated and the gills were perfused with sea water at low pressure while an irrigation system kept the fish at sea water temperature. For the next 3 h temperature measurements were made with thermistors. These probes were inserted into the muscle at a number of locations and isotherms were drawn by hand. Musyl et al.(Musyl et al., 2003) captured 87 bigeye tuna by conventional long-line gear and by rod and reel near the main Hawaiian Islands from April 1998 to April 2001. On board Wildlife Computers (Redmond, WA, USA) archival tags were implanted through a 1–1.5 cm incision into the peritoneal cavity and placed approximately as illustrated by the black circle in Fig. 1B. The tags sampled data on four channels every 60 s: depth,light, internal and external temperature. Post-operationally the fish were released and the data were downloaded from the tags from recaptured tuna.


To analyze the cross-sectional isotherms and the time-dependent temperature changes, we developed a 2D heat flux model of a cross-section of a bigeye tuna. The model included the following assumptions.

  1. The thermal properties of tissues/organs (muscle, bone, gut, etc.) differ but are identical everywhere within any tissue/organ.

  2. The given cross-section is representative for the entire tuna [i.e. the tuna is `infinite' in length (the z-coordinate) and does not alter shape].

  3. There is a non-directional, uniform convective heat transfer to the environment by the blood. This heat transfer is modified by the presence of vascular heat exchangers.

  4. There is no significant thermal boundary layer and the skin temperature is uniform and everywhere equal to the ambient water temperature.

  5. The heat flux of the cross-section is symmetrical across the vertical axis of the tuna (i.e. the net heat flux across this symmetry line is zero).

  6. All energy liberated in metabolism is eventually converted into heat.

Heat dissipation through the tuna (and thus through the cross-section)follows Fourier's law of heat conduction. This results in a second-order partial differential equation:
\[\ \frac{{\partial}T_{\mathrm{b}}}{{\partial}t}={\nabla}{\cdot}({\kappa}{\nabla}T_{\mathrm{b}})={\kappa}\left(\frac{{\partial}^{2}T_{\mathrm{b}}}{{\partial}x^{2}}+\frac{{\partial}^{2}T_{\mathrm{b}}}{{\partial}y^{2}}\right).\]
The last term in Eqn 1 holds for the 2D homogeneous case with constant thermal diffusivity. Inclusion of the metabolic heat produced within the body of the tuna and the convective loss by the blood, which is cooled in the gills, yields the following equation:
\[\ \frac{{\partial}T_{\mathrm{b}}}{{\partial}t}={\nabla}{\cdot}({\kappa}{\nabla}T_{\mathrm{b}})+{\dot{T}}_{\mathrm{p}}-{\varphi}(T_{\mathrm{b}}-T_{\mathrm{w}})(1-{\varepsilon}),\]
where p is temperature production, i.e. the rate of change in temperature due to heat production(°C s–1, see Appendix), Tb is body temperature at a given coordinate (°C), Tw is water temperature (°C), κ=K/(ρcp) is thermal diffusivity (cm2 s–1), K is thermal conductivity (J s–1 cm–1°C–1), ρ is density (g cm–3), cp is specific heat capacity (J g–1°C–1), β=ρcp is the heat capacitance coefficient (J cm–3 °C–1),ϵ is heat exchanger efficiency (dimensionless),φ=bbloodmuscle)is the convective cooling rate constant (s–1), b is specific blood flow(
), t is time, ▿ is the gradient operator, ▿· is the divergence operator. A derivation of the equation and a more detailed account of the meaning of the operators and the different terms are given in the Appendix.

Because second-order non-linear partial differential equations are virtually impossible to solve analytically, this was accomplished vianumerical analysis by the methods of finite elements using the commercial software package FlexPDE® (PDE Solutions, Inc., CA, USA).

The steady-state model

At a sea water temperature of 20°C (Tw), the 70 kg bigeye tuna presented in Carey and Teal(Carey and Teal, 1966) had a body temperature (Tb) of 32°C in the center of the red muscle. The highest Tb was found deep in the lateral muscle behind the pectoral fin, which is the thickest part of the fish. We developed a model to predict the isotherms of the cross-section of this fish.

We estimated body depth (the vertical distance, i.e. height of the fish in a transection at the thickest part of the body) relative to fork length and body mass by applying a geometric model of tuna body shape developed by Magnuson and Weininger (Magnuson and Weininger, 1978). We calculated fork length to be 144 cm (based on the relationship between body mass (M) and fork length (L)in bigeye tuna: M=0.00961L3.18) and body depth to be 42.8 cm. Once the scale was in place, the construction of the shape of the cross-section was accomplished using a spline function. To account for all thermally conductive regions of the body, the respective domains were implemented in the model (Fig. 1) and boundary conditions were specified. Each region was assigned a specific thermal property (Table 1). In a highly vascularized tissue such as the red muscle the natural convective fluctuations in the xy-coordinate direction due to blood flow could be expected to increase the effective thermal diffusivity,κ. However, experiments on tissue thermal conductivity performed on rats indicated that the contribution of blood flow to thermal diffusivity is negligible (Song et al., 1999;Cheng et al., 2002).

The volume-specific heat production in muscle was determined from the maximum oxygen consumption(O2max), 2700 mg O2 kg–1 h–1(Bushnell and Brill, 1991; Dewar and Graham, 1994),assuming all the energy extracted from the oxygen is eventually converted into metabolic heat (i.e. p=O2×calorific coefficient of O2). The maximum metabolic rate was selected in order to assess the maximal heat production, thus setting the lower limit for the estimated heat exchanger efficiencies. The volume-specific heat production in the red and white muscle was calculated from the average heat production and the relative area contributions of the two types of muscle by assuming their ratio was equal to the ratio of capillarity in the respective tissues(5.3≈2880 mm–2/543 mm–2) (from Mathieu-Costello et al.,1992). The specific blood flow in the convective cooling rate constant was calculated from:
\[\ {\dot{{\nu}}}_{\mathrm{b}}=\frac{{\dot{V}}_{\mathrm{O}_{2}}}{{\rho}{\Delta}S{\times}O_{2}{\ }\mathrm{carrying\ capacity}}.\]
The oxygen-carrying capacity was set at 0.16 ml O2
(Malte et al., 2007). The saturation difference between oxygenated and deoxygenated blood (arterial and venous blood) was set at ΔS=0.5. When solving Eqn 2, FlexPDE® utilized the finite element method of separating the domain in question into mesh grids of lengths that are finitely minute. This resulted in a linear approximation of the partial differential equation until a pre-selected margin of error was reached (Fig. 2).
To solve the heat transfer problem at steady state, Eqn 2 must equal zero. That is:
\[\ {\nabla}{\cdot}({\kappa}{\nabla}T_{\mathrm{b}})+{\dot{T}}_{\mathrm{p}}-{\varphi}(T_{\mathrm{b}}-T_{\mathrm{w}})(1-{\varepsilon})=\frac{{\partial}T_{\mathrm{b}}}{{\partial}t}=0.\]
Because of symmetry, only half the cross-section of the tuna was modeled by applying a boundary condition of zero heat flux across the line of symmetry(i.e. ∂Tb/∂n=0, where n is the normal to the surface). On the rest of the boundary the condition Tb=Tw was applied. By definition, the term in Eqn 4 will only approach zero at a given accuracy and an error limit was set at 1×10–6 for all calculations.

To simulate the observed isotherms obtained by Carey and Teal(Carey and Teal, 1966) at 20°C Tw, a reasonable heat exchanger efficiency value had to be estimated for both the red and white muscle. The estimated value was considered satisfactory if the modeled isotherms closely resembled the experimentally determined isotherms.

The time-dependent model

Bigeye tuna descend at dawn to depths of 300–500 m and make regular upward excursions into the warmer surface layer, then ascend at dusk to remain at the surface layer until the following dawn(Musyl et al., 2003). The upward excursions are made in order to warm the body rapidly. The fish then descends while retaining the newly gained heat. Fish 241 (from Musyl et al., 2003) was selected because it clearly demonstrated this vertical movement pattern characteristic of bigeye tuna (Fig. 3).

The non-steady-state model was based on a modification of the steady-state model. The time dependency relied on the fact that Tw (and therefore also Tb) varies in time as the tuna descends and ascends, thus altering the model into a non-steady-state problem. The processed external (Tw) temperature data from the archival tag of Fish 241 were imported into the model and these were used in Eqn 2 so that Tw=Tw(t) during its solution. The shape of the 2D cross-section of the bigeye tuna was set as identical to the steady-state model, but the dimensions were resized to fit the physical proportions of Fish 241. The calculations on the new dimensions were based on the formerly presented geometric model. Fish 241 had a fork length of 79 cm and a body mass of 10.3 kg, yielding a modeled body depth and body mass of 23.5 cm and 10.4 kg, respectively. The model was designed to allow for scaling with the placement of the internal organs preserved. A sensor region was introduced to the model to simulate the actual sensor position of the archival tag. As in Musyl et al. (Musyl et al.,2003) it was placed in the visceral cavity (Figs 1 and 2) and given an equivalent dimension.

An imitated shunt mechanism was implemented to allow for active thermoregulation. The shunt mechanism was assumed to direct the blood flow either: (1) through the countercurrent heat exchangers (hence conserving metabolically produced heat) or (2) around the countercurrent heat exchangers(possibly by passing the blood through the dorsal aorta and post-cardinal vein). This is a plausible route of active thermoregulation in the bigeye tuna as suggested by Malte et al. (Malte et al., 2007). We chose to use only the simplest possible mode of thermoregulation: either all the blood flow is through the exchanger or all the blood is bypassing the exchanger. This was achieved in the model by letting the efficiency (ϵ) of the heat exchanger vary according to:
\[\ {\varepsilon}={\varepsilon}_{\mathrm{fit}}(1-f),\]
where ϵfit is the heat exchanger efficiency fitted from the steady-state model and f is the fraction of blood bypassing the rete. The shunt fraction is then varied according to:
\[\ f=\begin{array}{ll}f_{\mathrm{high}}&\mathrm{for}{\ }T_{\mathrm{w}}-T_{\mathrm{b}}{>}0\\f_{\mathrm{low}}&\mathrm{for}{\ }T_{\mathrm{w}}-T_{\mathrm{b}}{<}0\end{array}.\]
Thus ϵ in Eqn 5 attains a high value whenever Tw<Tb, and a low value whenever Tw>Tb. The oxygen consumption was assumed to have either of two values – a`resting' and an active value (i.e. a two-state diurnal metabolism). When in nocturnal swimming mode, the bigeye tuna display a reduced vertical swimming activity, therefore referred to here as `resting'. They swim in the horizontal plan, possibly to catch prey, while maintaining the minimum speed required to ram ventilate and avoid sinking. At dawn, the fish enters an active swimming mode, descending to greater depths and then either feeding or swimming to the surface and down again, increasing its oxygen consumption. Calculations of the vertical swimming speed based on the depth data exhibited a clear consistency with the assumption of low cruising mode at night with low swimming speeds,and an active swimming mode during the day with high swimming speeds(Fig. 3). Hence the shift in oxygen consumption was modeled based on the day/night cycle, simply as:
\[{\dot{V}}_{\mathrm{O}_{2}}(x)=\begin{array}{ll}{\dot{V}}_{\mathrm{O}_{2}\mathrm{max}}&\mathrm{for}{\ }x=1\\{\dot{V}}_{\mathrm{O}_{2}\mathrm{rest}}&\mathrm{for}{\ }x=0,\end{array}\]
where x is an array in which daytime is given the number 1 and 0 is for night-time. O2rest was set at a routine metabolic rate of 540 mg O2 kg–1h–1 (Freund,1999). Thermal properties, heat production p and the fitted heat-conserving efficiency (ϵfit) for red and white muscle from the steady-state model were applied. The model was considered satisfactory if it predicted values of Tb at the sensor region that intimately followed the Tb values recorded by the sensor tag and the adjusted parameter was now the shunt fraction(Eqn 5).

Fig. 4 illustrates cross-sections of the recorded isotherms(Fig. 4A) and the isotherms predicted by the model (Fig. 4B). It is obvious that the observed and predicted results match and that the model simulates the 2D heat flux in the fish quite well. The highest temperatures are in the heat-generating red muscle from which there is a heat flux down the thermal gradient into the adjacent white muscle. A fraction of the heat flux will pass through the gut and finally reach the surface of the bigeye tuna in contact with 20°C sea water. The best match,as shown in the figure, was obtained when the efficiency of the heat exchangers supplying blood to the red and white muscle was 99% and 96%,respectively.

To better display the precision of our estimates, Fig. 5 shows the temperature profile along the secant line through the isotherms of the cross-section illustrated in Fig. 4B. This specific cross-section includes all observed isotherms traversing the point of maximal temperature in the red muscle (Fig. 5C). Fig. 5Adisplays the smoothly fitting curvature of the predicted body temperature profile and the isotherms recorded by Carey and Teal(Carey and Teal, 1966). This close fit was only achieved using the above-mentioned heat exchanger efficiencies. The deviations of the modeled Tb profile are not considerably different from the observed isotherms, when disregarding the outlier at the 26°C isotherm. Fig. 5B shows the significance of a high degree of heat conservation in the white muscle. It is clear that only at ϵwhite=96% is it possible to recreate the observed isotherms (ϵred being fixed at 99%), thus stressing the importance of heat conservation in the white muscle. The impact of the high efficiency rete in the blood supply to the white muscle can be seen when comparing the 90% curve (blue) and 96% curve(red). The integrated area under the 96% curve is almost twice as large(235/125=1.88) as the area under the 90% curve.

To further explore the consequence of an absence of a rete in the blood supply to the white muscle (ϵwhite=0%), we plotted the temperature distribution in the 2D cross-section(Fig. 6). The temperature distribution is very different from that in Fig. 4B. Any heat produced in the white muscle will rapidly be lost via the convective blood flow resulting in `lumps' of isotherms at the boundary of the red muscle. Even though the efficiency of the rete in the red muscle is maintained at 99%, the warmest point is only 26°C as opposed to 32°C in Fig. 4. We thus conclude that the high temperatures and isothermal distribution observed by Carey and Teal(Carey and Teal, 1966) require the presence of an effective countercurrent heat exchanger in the vascular system supplying blood to white muscle.

To quantify the loss of heat across the surface of the fish, the heat flux vectors (Fig. 7A) at the periphery were integrated (Fig. 7B). This yielded a surface heat loss of 0.622 J s–1 cm–1 at a total heat production of 1.886 J s–1 cm–1 (i.e. surface heat loss accounts for approximately 33% of total heat production). Despite the relatively high efficiencies of the vascular countercurrent heat exchanger(ϵred=99%, ϵwhite=96%), approximately two-thirds of the heat produced is still lost via the gills. Ifϵ red and ϵwhite were 100%, all the heat produced would have to be lost across the surface. When testing this notion in the model, we found a 99.6% surface heat loss, thus confirming the suitability of our methods.

The mass to heat loss relationship was modeled to determine the significance of body size. This was accomplished via the scalability of the modeled cross-section of the bigeye tuna. The relationship between mass and surface heat loss:total heat production ratio is presented in Fig. 8 (blue curve). The hyperbolic shape of the curve is expected because of the surface to volume ratio (i.e. small fish have larger relative surface areas than larger fish). At a mass of approximately 7 kg, 50% of the heat production is lost via convection through the blood and 50% via conduction through the skin. In order to investigate the implication of an increased body size on the temperature in the warmest area of the body, the relationship between maximal Tb in the red muscle and mass was established (red curve in Fig. 8). The curve in the logarithmic plot is sigmoid and appears to converge to a limit [data not shown at an unrealistically high body mass:bigeye tuna 170.6 kg, bluefin tuna (Thunnus thynnus) 678.6 kg; http://www.bigmarinefish.com]. The curve fits well with the mass and temperature of the red muscle of the steady-state model. Additionally, the curve fits with isotherms in cross-sections of smaller fish throughout the body of the tuna examined by Carey and Teal (Carey and Teal,1966) (data not shown). A 1 kg fish would have a sustained Tx (maximal temperature elevation compared with the ambient water) of approximately 3.5°C, and a sustained Tx of 10°C required a mass of approximately 20 kg.

We used the archival tag data for a selected daytime period of Fish 241 of Musyl et al. (Musyl et al.,2003) to test the time-dependent model. The specific day was chosen based on the criterion that it was representative of the bigeye tuna behavioral swimming pattern (i.e. regular excursions in between the uniform-temperature surface layer and depths at approximately 350 m such as day 4 in Fig. 3). A typical example of water temperature, recorded body temperature, and modeled body temperature during the daytime for Fish 241 is displayed in Fig. 9A. Twrapidly decreased as the fish descended, while Tbdecreased much more slowly. For each excursion to the uniform-temperature surface layer, there was a rapid increase in Tbaccompanied by a slow decrease when returning to depth. This difference in heating and cooling rate considerably magnified the average difference between Tw and Tb. The average Tb is approximately 19.3°C during the regular excursions. However, the range of Tb during the regular excursions (approximately 17–22°C) is far more interesting. The tuna seems tolerant of a reduction in Tb to only 17°C. When Tb reached this value the tuna initiated the rapid ascent to warmer surface waters. Likewise, when Tb reached about 22°C, the fish commenced its descent. The modeled Tbwhich, for this purpose, was the average temperature of the sensor region in the geometry (Fig. 1), follows the recorded Tb remarkably well and only during the last three excursions does a minor, but noticeable, discrepancy develop. This close fit was achieved by letting the shunting fraction (f) vary simply between 0 and 1 (i.e. all the blood either passed through the heat exchangers or completely bypassed it; Eqn 5). Since Fish 241 used for the time-dependent model was somewhat smaller than the tuna used for the steady-state model, the efficiency of the heat exchanger supplying blood to the red muscle had to be increased from 99%to 99.5% to reproduce measured Tb values. The precision of the fit can be better appreciated in Fig. 9B, which shows an expanded segment of Fig. 9A. This figure reveals that rapid modification of the heating and cooling rates (physiological thermoregulation) can be achieved by shunting blood through or around the heat exchanger resulting in a significant degree of control over Tb. When preventing physiological thermoregulation by holding the parameters f and ϵ constant, the modeled Tb did not in any way resemble the measured Tb values (data not shown). This re-emphasizes similar findings previously obtained by Holland and Sibert(Holland and Sibert, 1993) and Malte et al. (Malte et al.,2007). We tested the effects of pure thermal inertia in the time-dependent model by scaling up the geometry, reducing the efficiency of the rete, and lowering the oxygen consumption; thus simulating a large ectothermic teleost. The result was a curve radically different from the simulations derived from modeling a thermoregulating bigeye tuna. In the model simulating a large ectothermic teleost, Tb follows Tw when the fish resides in the surface layer. However,when the fish descends into colder water, Tb follows Tw but with a noticeable delay. Similarly when the fish ascends there is a delay in Tb since the thermal inertia works both ways. However, even with frequent regular vertical excursions, the virtual ectotherm maintained a Tx of only 2°C (data not shown).

Steady-state heat exchange

It has been assumed that white muscle is warmed via conduction from the deep red muscle (Korsmeyer and Dewar, 2001). Based on our results, however, this appears not to be entirely correct. When we set the parameters in our model to assume that there is no rete in the vascular system supplying blood to the white muscle(i.e. zero rete efficiency and that the white muscle is warmed solely by heat conducted from the red muscle), the model does not predict the isotherms reported by Carey and Teal (Carey and Teal, 1966) (Fig. 6). Moreover, not even when setting ϵred=99.999%could we reproduce the recorded temperatures when ϵwhite=0. Only when multiplying the maximal metabolic rate by at least a factor 10 is it possible to reach red muscle temperatures 12°C above Tw.

Rather, our model shows that: the presence of a functional rete in the blood supply to the white muscle is necessary to achieve realistic model outputs, less than 33% of the heat produced by the red muscle is actually conducted to the white muscle, and the white muscle is warm primarily because of the heat that it itself produces, as well as the heat it receives from the red muscle. Moreover, the white muscle aids temperature elevation in the red muscle by exerting an insulating effect (i.e. lowering the thermal gradient between the red and white muscle) (Figs 4 and 6). The heat-conserving capacity of the white muscle therefore not only serves to elevate the metabolic activity of this tissue but also provides the potential for the red muscle to reach the high temperatures that Carey and Teal(Carey and Teal, 1966)recorded. The parameter ΔS (oxygen saturation difference between arterial and venous blood) was set at a value of 0.5. It could be argued that the ΔS in the muscle is higher during elevated swimming speeds. However, setting ΔS=1 yielded only a minor reduction in heat exchanger efficiency (ϵred=97% andϵ white=94%) compared with when ΔS=0.5.

The vascular countercurrent heat exchangers in the red muscle are indeed well developed with 4-to-1 modules arranged in tight compaction(Stevens et al., 1974). However, Carey and Teal (Carey and Teal,1966) described the rete in the white muscle as vascular bands near the skin surface which ultimately split up in triads in the deep muscle. Considering the clearly simpler construction of the rete in the blood supply to the white muscle, the ϵwhite=96% implied by the results of our calculations seems surprising. Preliminary simulations of heat fluxes in different rete arrangements, however, indicate that this is indeed possible(H.M., in preparation). Nevertheless, the vascular countercurrent heat exchanger in the red muscle still conserves heat four times better than the simpler triad heat exchanger in the white muscle (the factor 1–ϵ in Eqn 2).

Because of the larger mass of the white muscle, total heat production within this tissue is approximately 2.3 times larger than that in the red muscle. The ratio of red and white muscle mass to total body mass is 6–9% and 45–55%, respectively(Graham et al., 1983; Freund, 1999). Integrating the muscle mass in relation to the total body mass in the specific cross-section in our model yields a relative red and white muscle mass of 6.6% and 82.2%,respectively. The modeled cross-section is located just behind the pectoral fin, the thickest part of the fish, whereas the 45–55% is an expression of the total white muscle mass in all three dimensions. This explains the discrepancy in the fraction of white muscle mass.

The finding of the white muscle as a significant contributor in heat conservation is also relevant to the evolution of endothermy in fish. Our results suggest that the rete in the blood supply to the white muscle most likely evolved either simultaneously with, or before, the rete in the blood supply to the red muscle and the internalization of the red muscle. If it had not, our results clearly imply that internalization of the red muscle alone would not allow significant sustained temperature elevation of the red muscle without the presence of a rete in the blood supply to the white muscle.

The time-dependent setting

In the time-dependent model, the assumed blood flow patterns were either that all the blood went through the rete to conserve metabolically produced heat (when Tb>Tw) or all the blood bypassed the rete (when Tw>Tb). Likewise, the oxygen consumption (leading to heat production) was modeled in a simple two-state fashion: constant at maximal metabolic rate during the daytime, and constant at routine metabolic rate at night-time. In spite of this simple all-or-none arrangement, our model reproduced the recorded body temperature data very well. The all-or-none shunting mechanism was previously suggested by Malte et al. (Malte et al.,2007). We may seem less justified in simplifying the modeled diurnal metabolism to a two-state oxygen consumption rate. However, analysis of the data of Fish 241 revealed similar velocities of ascents and descents with a gradual increase in swimming speed towards mid-point depth followed by a gradual decrease (Malte et al.,2007). Holland and Sibert(Holland and Sibert, 1993)made similar assumptions. Moreover, there was only a minor difference [<8%(Malte et al., 2007)] between cooling and heating rates in a lumped parameter model when oxygen consumption was held constant, or when oxygen consumption was varied with vertical swimming speed during vertical excursions. Even though a simple all-or-none mode of regulation explains the present data well a more finely graded regulation can certainly not be excluded. Testing and comparing different modes of regulation would require more detailed knowledge of activity (not only vertical swimming velocity) and temperature data of better spatial and temporal resolution than are available today.

The cooling and heating rate constants found in the lumped parameter model during the vertical excursions were in the range 0.02–0.04 min–1 and 0.2–0.6 min–1, respectively. When we estimated these from the modeled body temperature, we arrived at similar rates (data not shown). Therefore, the present model shows that these cooling and heating rate constants can indeed be the result of an all-or-none shunting mechanism (e.g. diverting blood from the lateral heat exchangers to the dorsal aorta/post-cardinal vein and vice versa).

The time-dependent model predicted the recorded Tb very well for various fish on different days. However, in a few fish the metabolic heat production in the white muscle had to be multiplied by a factor of two to yield results that resembled the observed data. This may be because of significant anaerobic metabolism under certain conditions. Tuna white muscle is known to have a large anaerobic capacity with high glycogen stores and extraordinarily high lactate dehydrogenase activity(Hulbert et al., 1979; Arthur et al., 1992; Dickson, 1995). The model only takes aerobic heat production into account and ignores heat production associated with anaerobic glycolysis thus missing the actual scope of white muscle heat production. Increasing white muscle heat production would imply a corresponding decrease in the necessary efficiency of the countercurrent heat exchanger in the blood supply to the white muscle. However, this would only be temporary and, even then, modeling of all realistic scenarios still demonstrated the necessity of a white muscle rete efficiency of around 95%.

Effects of body mass

Scaling the model tuna makes it apparent that body mass is a major contributor to elevated body temperature, thus supporting the similar findings of Malte and colleagues (Malte et al.,2007). Assuming Tw=20°C, a constant volume-specific oxygen consumption rate and high rete efficiencies as estimated from the present calculations, our model predicts that very large tuna would overheat. Thus, in this setting, a 400 kg tuna would reach a core body temperature of approximately 40°C. This, like in all endotherms,stresses the importance of a mass-specific metabolic rate that decreases with size. At this point it is difficult to draw specific conclusions with regard to the value of the exponent of the metabolic rate–body mass relationship of tuna because of the limited size range for which there are data (Korsmeyer and Dewar,2001). It is clear from Fig. 8, however, that fractional surface heat loss decreases as body size increases and the significance of heat exchanger efficiency in controlling heat loss increases. Thus, a size-dependent heat loss viathis route could be achieved by a size-dependent diameter of the arterioles and venules comprising the heat exchangers. There may be an indication of this since Stevens and colleagues (Stevens et al., 1974) found the arteriole and venule diameters to be around 40 and 80 μm in the (central) heat exchanger of a 1.9 kg skipjack tuna whereas Carey and Teal (Carey and Teal,1966) reported values of around 100 and 200 μm in the similarly arranged (lateral) heat exchanger of a 70 kg bigeye tuna. However, such anatomical data are too scarce to allow definitive conclusions to be drawn. Apart from increasing the possible steady-state Tx, a larger body mass also will slow cooling rates and thus decrease the frequency of vertical excursions.

In conclusion, we have shown that an effective vascular countercurrent heat exchanger in the blood supply of both the red and white muscle is a prerequisite for endothermy in tuna. Furthermore, physiological thermoregulation, as expressed by differential cooling and heating rates, can be fully accounted for by shunting of blood either through or around vascular heat exchangers in an all-or-none fashion. Moreover, by adjusting the heat exchanger efficiencies obtained from steady-state modeling of a large bigeye tuna, we were able to predict the temperature variations recorded in a much smaller tuna performing daily vertical excursions. This suggests that the estimated heat exchanger efficiencies are fairly conservative and that the higher body temperatures in larger individuals are primarily due to a lower fractional surface heat loss and a greater thermal inertia.

Suggestions for future research

This study points to the necessity for studies on heat conservation in white muscle. As a first step, it will be necessary to obtain quantitative anatomical data (i.e. vessel diameter and length as well as interspacing of heat exchanger units) in vascular bands and triads serving the white muscle. This will make it possible to model the maximum possible efficiency of most idealized rete arrangements.

Model development

Eqn A1 is based on a heat balance on an infinitesimal element. That is, the heat generated within a unit volume of tissue, the heat flux (i.e. conduction) in all three coordinate directions, and the change in energy per unit time must balance:
\[\ q_{x}+q_{y}+q_{z}+q_{\mathrm{gen}}=q_{x+\mathrm{d}x}+q_{y+\mathrm{d}y}+q_{z+\mathrm{d}z}+\frac{\mathrm{d}E}{\mathrm{d}t},\]
where qx, qy and qz are the heat fluxes into the unit volume in the x-, y- and z-directions (respectively), qgen is the heat generated per unit time within the unit volume, qx+dx,qy+dy and qz+dz are the heat fluxes out of the unit volume in the x-, y- and z-directions (respectively) and dE/dt is the rate of change of energy content in the unit volume.
These energy quantities are given as follows:
\[\ \begin{array}{ll}&q_{x}=-K{\ }\mathrm{d}y{\ }\mathrm{d}z\frac{{\partial}T}{{\partial}x},\\&q_{x+\mathrm{d}x}=-\left[K\frac{{\partial}T}{{\partial}x}+\frac{{\partial}}{{\partial}x}\left(K\frac{{\partial}T}{{\partial}x}\right)\mathrm{d}x\right]\mathrm{d}y{\ }\mathrm{d}z,\end{array}\]
\[\ \begin{array}{ll}&q_{y}=-K{\ }\mathrm{d}x{\ }\mathrm{d}z\frac{{\partial}T}{{\partial}y},\\&q_{y+\mathrm{d}y}=-\left[K\frac{{\partial}T}{{\partial}y}+\frac{{\partial}}{{\partial}y}\left(K\frac{{\partial}T}{{\partial}y}\right)\mathrm{d}y\right]\mathrm{d}x{\ }\mathrm{d}z,\end{array}\]
\[\ \begin{array}{ll}&q_{z}=-K{\ }\mathrm{d}x{\ }\mathrm{d}y\frac{{\partial}T}{{\partial}z},\\&q_{z+\mathrm{d}z}=-\left[K\frac{{\partial}T}{{\partial}z}+\frac{{\partial}}{{\partial}z}\left(K\frac{{\partial}T}{{\partial}z}\right)\mathrm{d}z\right]\mathrm{d}x{\ }\mathrm{d}y,\end{array}\]
\[\ q_{\mathrm{gen}}={\dot{h}}_{\mathrm{p}}{\ }\mathrm{d}x{\ }\mathrm{d}y{\ }\mathrm{d}z,\]
\[\ \frac{\mathrm{d}E}{\mathrm{d}t}={\rho}c_{\mathrm{p}}{\ }\mathrm{d}x{\ }\mathrm{d}y{\ }\mathrm{d}z\frac{{\partial}T}{{\partial}t},\]
where K is the thermal conductivity of the material with the units J s–1 cm–1 °C–1, T is the temperature as a function of a given position and time, p is the volume-specific metabolic heat production in J s–1 cm–3,ρ is the density with the units g cm–3 and cp is the specific heat capacity with the units J g–1 °C–1.
Inserting Eqns A2, A3, A4, A5, A6 into Eqn A1 yields:
\[\ {\rho}c_{\mathrm{p}}\frac{{\partial}T}{{\partial}t}=\frac{{\partial}}{{\partial}x}\left(K\frac{{\partial}T}{{\partial}x}\right)+\frac{{\partial}}{{\partial}y}\left(K\frac{{\partial}T}{{\partial}y}\right)+\frac{{\partial}}{{\partial}z}\left(K\frac{{\partial}T}{{\partial}z}\right)+{\dot{h}}_{\mathrm{p}}.\]
With the assumption of uniformity in the z-direction, the differential term in the z-direction in Eqn A7 vanishes. We now introduce a non-directional, uniform, convective cooling term. This can be thought of as a uniform blood flow perpendicular to the xy-plane. This convective flow through the tissue is assumed to have a cooling power at any point proportional to the temperature difference between that point and the surrounding water (TbTw). Furthermore it is assumed to be modified by the presence of a heat exchanger with a factor (1–ϵ) where ϵ is the heat exchanger efficiency. This leads to:
\[\ {\rho}c_{\mathrm{p}}\frac{{\partial}T}{{\partial}t}=\frac{{\partial}}{{\partial}x}\left(K\frac{{\partial}T}{{\partial}x}\right)+\frac{{\partial}}{{\partial}y}\left(K\frac{{\partial}T}{{\partial}y}\right)+{\dot{h}}_{\mathrm{p}}+{\dot{{\nu}}}_{\mathrm{b}}{\beta}_{\mathrm{b}}(T_{\mathrm{w}}-T_{\mathrm{b}})(1-{\varepsilon}).\]
By dividing through with ρcp (for muscle), and introducing thermal diffusivityκ=K/(ρcp) and temperature production p=p/(ρcp)we get:
\[\ \frac{{\partial}T}{{\partial}t}=\frac{{\partial}}{{\partial}x}\left({\kappa}\frac{{\partial}T}{{\partial}x}\right)+\frac{{\partial}}{{\partial}y}\left({\kappa}\frac{{\partial}T}{{\partial}y}\right)+{\dot{T}}_{\mathrm{p}}+{\dot{{\nu}}}_{\mathrm{b}}\frac{{\beta}_{\mathrm{b}}}{{\beta}_{\mathrm{m}}}(T_{\mathrm{w}}-T_{\mathrm{b}})(1-{\varepsilon}).\]
Finally, introducing the gradient (▿) and the divergence (▿·)operators and defining a rate constant (dimension s–1) for convective coolingφ=bβbmwe arrive at:
\[\ \frac{{\partial}T}{{\partial}t}={\nabla}{\cdot}({\kappa}{\nabla}T)+{\dot{T}}_{\mathrm{p}}+{\varphi}\left(T_{\mathrm{b}}-T_{\mathrm{w}}\right)(1-{\varepsilon}).\]
This is the second-order partial differential equation at its final edition,which is a key element in the models. For constant thermal diffusivity, the first term on the right-hand side of Eqn A10 reduces to κ▿2Twhere▿2=▿·▿ is the Laplace operator. This is the case inside every homogeneous domain but not on its boundaries.

The last term (the convective term) is thus dependent on the difference between Tb and Tw, and the efficiency of the heat exchanger. Tb will be depend on x and y as well as on time [i.e. Tb=Tb(x,y,t)], whereas Tw may or may not depend on time but does not depend on x or y. The efficiency takes on values in the interval 0–1 (i.e. a large ϵ results in a reduced convective loss). In other words, the final partial differential equation(Eqn A10) is an expression for a 3D model modified into a 2D model taking into account the z-dimension as a uniform convective blood flow perpendicular to the xy-plane.


  • cp

    specific heat capacity (J g–1°C–1)

  • f

    shunting fraction

  • p

    volume-specific heat production (J s–1cm–3)

  • K

    thermal conductivity (J s–1 cm–1°C–1)

  • L

    fork length

  • M

    body mass

  • n

    normal to the surface

  • S

    O2 saturation

  • t


  • Tb

    body temperature at a given coordinate (°C)

  • p

    temperature production, i.e. the rate of change in temperature due to heat production (°C s–1)

  • Tw

    water temperature (°C)

  • Tx

    maximal temperature elevation compared with the ambient water(°C)

  • b

    specific blood flow(


  • O2max

    maximum oxygen consumption (mg O2 kg–1h–1)

  • O2rest

    resting oxygen consumption (mg O2 kg–1h–1)

  • gradient operator

  • ▿·

    divergence operator

  • ▿2

    Laplace operator

  • β

    heat capacitance coefficient (β=ρcp in J cm–3 °C–1)

  • ϵ

    heat exchanger efficiency (dimensionless)

  • κ

    thermal diffusivity [κ=K/(ρcp)in cm2 s–1]

  • ρ

    density (g cm–3)

  • φ

    convective cooling rate constant[(φ=bbloodmuscle)in s–1]

This work was supported by The Pelagic Fisheries Research Foundation (PFRP) and The Danish Natural Sciences Research Foundation (FNU).

Altringham, J. D. and Shadwick, R. E. (
). Swimming and muscle function. In
Tuna: Physiology, Ecology And Evolution, Fish Physiology
, vol.
(ed. B. A. Block and E. D. Stevens), pp.
-341. San Diego: Academic Press.
Arthur, P. G., West, T. G., Brill, R. W., Schulte, P. M. and Hochachka, P. W. (
). Recovery metabolism of skipjack tuna(Katsuwonus pelamis) white muscle: Rapid and parallel changes in lactate and phosphocreatine after exercise.
Can. J. Zool.
Blank, J. M., Morrissette, J. M., Landeira-Fernandez, A. M.,Blackwell, S. B., Williams, T. D. and Block, B. A.(
). In situ cardiac performance of Pacific bluefin tuna hearts in response to acute temperature change.
J. Exp. Biol.
Block, B. A., Finnerty, J. R., Stewart, A. F. R. and Kidd,J. (
). Evolution of endothermy in fish: mapping physiological traits on a molecular phylogeny.
Bonnick, S. L. (
Osteoporosis,The Hand Book, third edition
, 147 pp. Texas: Cooper Square Press.
Brill, R. W. and Bushnell, P. G. (
). The cardiovascular system of tunas. In
Tuna: Physiology, Ecology And Evolution, Fish Physiology
, vol.
(ed. B. A. Block and E. D. Stevens), pp.
-120. San Diego: Academic Press.
Bushnell, P. G. and Brill, R. W. (
). Responses of swimming skipjack (Katsuwonus pelamis) and yellowfin(Thunnus albacares) tunas to acute hypoxia and a model of their cardiorespiratory function.
Physiol. Zool.
Calttenburg, R., Cohen, J., Conner, S. and Coo, N.(
). Thermal properties of cancellous bone.
J. Biomed. Mater. Res.
Carey, F. G. (
). Fishes with warm bodies.
Sci. Am.
Carey, F. G. and Teal, J. M. (
). Heat conservation in tuna fish muscle,
Proc. Natl. Acad. Sci. USA
Carey, F. G., Teal, J. M., Kanwisher, J. W. and Lawson, K. D. (
). Warm bodied fish.
Am. Zool.
Carey, F. G., Teal, J. M., Kanwisher, J. W. and Stevens, E. D. (
). Bluefin tuna warm their viscera during digestion.
J. Exp. Biol.
Cheng, H. M. and Plewes, D. B. (
). Tissue thermal conductivity by magnetic resonance thermometry and focused ultrasound heating.
J. Magn. Reson. Imaging
Dewar, H. and Graham, J. B. (
). Studies of tropical tuna swimming performance in a large water tunnel, I. Energetics.
J. Exp. Biol.
Dewar, H., Graham, J. B. and Brill, R. W.(
). Studies of tropical tuna swimming performance in a large water tunnel, II. Thermoregulation.
J. Exp. Biol.
Dickson, K. A. (
). Unique adaptations of the metabolic biochemistry of tunas and billfishes for life in the pelagic environment.
Environ. Biol. Fishes
Freund, E. V. (
Comparisons of metabolic and cardiac performance in scombrid fishes: Insights into the evolution of endothermy
. PhD dissertation, Biological Sciences,Stanford University, Stanford, CA.
Graham, J. B. and Dickson, K. A. (
). The evolution of thunniform locomotion and heat conservation in scombrid fishes:new insights based on the morphology of Allothunnus fallai.
Zool. J. Linn. Soc. Lond.
Graham, J. B. and Dickson, K. A. (
). Anatomical and physiological specializations for endothermy. In
Tuna: Physiology, Ecology And Evolution, Fish Physiology
, vol.
(ed. B. A. Block and E. D. Stevens), pp.
-165. San Diego: Academic Press.
Graham, J. B., Koehrn, F. J. and Dickson, K. A.(
). Distribution and relative proportions of red muscle in scombrid fishes: Consequences of body size and relationships to locomotion and endothermy.
Can. J. Zool.
Gunn, J. and Block, B. A. (
). Advances in acoustic, archival and satellite tagging of tunas. In
Tuna:Physiology, Ecology And Evolution, Fish Physiology
, vol.
(ed. B. A. Block and E. D. Stevens), pp.
-224. San Diego: Academic Press.
Hensel, H. and Bock, K. D. (
). Durchblutung und Wärmeleitfährigkeit des menschlichen Muskels.
Phlügers Arch.
Holland, K. N. and Sibert, J. R. (
). Physiological thermoregulation in bigeye tuna, Thunnus Obesus.
Environ. Biol. Fishes
Holman, J. P. (
Heat Transfer
. Singapore: McGraw-Hill.
Hulbert, W. C., Guppy, M., Murphy, B. and Hochachka, P. W.(
). Metabolic sources of heat and power in tuna muscles. I. Muscle fine structure.
J. Exp. Biol.
Katz, S L., Syme, D. A. and Shadwick, R. E.(
). High-speed swimming, Enhanced power in yellowfin tuna.
Korsmeyer, K. E., and Dewar, H. (
). Tuna metabolism and energetics. In
Tuna: Physiology, Ecology and Evolution, Fish Physiology
, vol.
(ed. B. A. Block and E. D. Stevens), pp.
-78. San Diego:Academic Press.
Lowe, T. E., Brill, R. W. and Cousins, K. L.(
). Blood oxygen-binding characteristics of bigeye tuna(Thunnus obesus), a high-energy-demand teleost that is tolerant to low ambient oxygen.
Mar. Biol.
Magnuson, J. J. and Weininger, D. (
). Estimation of minimum sustained speed and associated body drag of scombrids. In
The Physiological Ecology of Tunas
, Academic Press,New York.
Malte, H., Larsen, C., Musyl, M. and Brill, R.(
). Differential heating and cooling rates in bigeye tuna(Thunnus obesus, Lowe): a model of non-steady-state heat exchange.
J. Exp. Biol.
Mathieu-Costello, O., Agey, P. J., Logemann, R. B., Brill, R. W. and Hochachka, P. W. (
). Capillary-fiber geometrical relationships in tuna red muscle.
Can. J. Zool.
Musyl, M. K. Brill, R. W., Boggs, C. H., Curran, D. S., Kazama,K. K. and Ski, M. P. (
). Vertical movements of bigeye tuna (Thunnus obesus) associated with islands, buoys, and seamounts near the main Hawaiian Islands from archival tagging data.
Fish. Oceanogr.
Ngadi, M. O. and Ikediala, J. N. (
). Heat transfer properties of chicken-drum muscle.
J. Sci. Food Agric.
Sepulveda, C. A., Dickson, K. A. and Graham, J. B.(
). Swimming performance studies on the eastern Pacific bonito (Sarda chiliensis), a close relative of the tunas (Family scombridae). I. Energetics.
J. Exp. Biol.
Song, J., Xu, L. X., Lemons, D. E. and Weinbaum, S.(
). Microvascular thermal equilibration in rat spinotrapezius muscle.
Ann. Biomed. Eng.
Stevens, E. D. and McLeese, J. M. (
). Why bluefin tuna have warm tummies: temperature effect on trypsin and chymotrypsin.
Am. J. Physiol.
Stevens, E. D., Lam, H. M. and Kendall, J.(
). Vascular anatomy of the counter-current heat exchanger of Skipjack Tuna.
J. Exp. Biol.
Urbanchek, M. G., Picken, E. B., Kalliainen, L. K. and Kuzon,W., Jr (
). Specific force deficit in skeletal muscles of old rats in partially explained by the existence of denervated muscle fibers.
J. Gerontol. A Biol. Sci. Med. Sci.