PO2 gradiants
import numpy as np
import pandas as pd
import statsmodels.api
import scipy.stats
import bokeh.io
import bokeh.plotting
#local .py file for some plotting functions and non-parametric bootstrapping utils
import plotting_utils
import numba
bokeh.io.output_notebook()
We can now read in the data into a dataframe for analysis.
df = pd.read_csv("./20190322_supp_table_2.csv")
We will add some calculations/generate a species averaged version of the dataframe.
df['species_underscore'] = [spec.replace(" ", "_") for spec in df['species']]
df_averages = df.groupby(['species', 'species_underscore', 'spiracle'], as_index=False).aggregate(np.average)
df_averages['subfamily'] = df.groupby(['species', 'species_underscore', 'spiracle'], as_index=False).aggregate(max)['subfamily']
species_per_subfam=df_averages.groupby(['subfamily', 'spiracle'], as_index=False).count().groupby('subfamily').aggregate(max).reset_index()[['subfamily', 'species']]
species_per_subfam.columns = ('subfamily', 'subfam_count')
df_averages = df_averages.merge(species_per_subfam, on='subfamily')
df_averages['log area (mm^2)'] = np.log10(df_averages['area (mm^2)'])
df_averages['log dist'] = np.log10(df_averages['depth (mm)'])
df_averages['log mass (g)'] = np.log10(df_averages['mass (g)'])
df_averages['log area/dist'] = np.log10(df_averages['area (mm^2)']/df_averages['depth (mm)'])
df_averages['log area^2/dist'] = np.log10(df_averages['area (mm^2)']**2/df_averages['depth (mm)'])
df['log area (mm^2)'] = np.log10(df['area (mm^2)'])
df['log dist'] = np.log10(df['depth (mm)'])
df['log mass (g)'] = np.log10(df['mass (g)'])
df['log area/dist'] = np.log10(df['area (mm^2)']/df['depth (mm)'])
df['log area^2/dist'] = np.log10(df['area (mm^2)']**2/df['depth (mm)'])
df_averages.head()
Using these morphological properties of the spiracles that we measured, we can calculate real world meaningful gas exchange properties with some unit conversions/physical gas properties. To calculate the diffusive capacity of the spiracles, we use the following relation:
\begin{align} G_{\mathrm{diff}} = \frac{\text{area}}{\text{depth}} * D * \beta \end{align}with $D$ (the diffusivity constant for O2 in air) $= 0.178 \text{ cm}^2$ and $\beta$ (the capacitance coefficient for oxygen in air) $= 404 \frac{\text{ nmol}}{\text{cm}^{3} \text{kPa}}$. With that value, you can calculate the oxygen consumption rate in $\frac{\text{nmol}}{\text{sec}}$ by multiplying the oxygen partial pressure gradient $\Delta \text{PO}_2$ across the spiracle in kPa by $G_{\text{diff}}$, that is
\begin{align} \text{Oxygen consumption rate }\left(\frac{\text{nmol}}{\text{sec}}\right) = \Delta \text{PO}_2 * G_{\mathrm{diff}} \end{align}To calculate advective gas transport capacities, we use $G_{\text{adv}}$, which comes from Poiseulle’s law and is
\begin{align} G_{\mathrm{adv}} = \frac{\text{area}^2}{8 * \pi * \text{dynamic viscosity} * \text{depth}} \end{align}where the dynamic viscosity of air is $1.86 * 10^{-8} \text{kPa sec}$. With these relations in hand, we can calculate both $G_{\text{diff}}$ and $G_{\text{adv}}$ for the beetle spiracles, both individually, and also by summing the individual capacities for all spiracles in the beetle.
df_test = df_averages.copy()
df_test['area/depth'] = (df_averages['area (mm^2)']/df_averages['depth (mm)'])/10
df_test['area^2/depth'] = ((df_averages['area (mm^2)']**2/df_averages['depth (mm)']))/(1000*1000*1000) #convert to cubic meters
df_test['g_diff'] = df_test['area/depth']*0.178*404
df_test['g_adv'] = df_test['area^2/depth']*(1/(np.pi*8*1.86*10**(-8)))
df_summed = df_test.groupby('species').median().reset_index()[['species', 'log mass (g)', 'g_diff', 'g_adv']]
df_summed['g_diff'] = df_test.groupby('species').sum().reset_index()['g_diff']
df_summed['g_adv'] = df_test.groupby('species').sum().reset_index()['g_adv']
df_summed['area'] = df_test.groupby('species').sum().reset_index()['area (mm^2)']
df_summed['average depth'] = df_test.groupby('species').mean().reset_index()['depth (mm)']
plots = []
resample_size = 10_000
lw = 2
cs = 12
m = df_summed['log mass (g)']
g_diff = np.log10(df_summed['g_diff']*2) #double the number for g_diff since we only measured one side of the animal, they have two of each spiracle
slope, intercept = np.polyfit(m, g_diff, deg=1)
x = np.array([m.min(), m.max()])
y = slope * x + intercept
p = bokeh.plotting.figure(plot_height=300, plot_width=400,
x_range=(x[0]-0.1, x[1]+0.1), y_range=(g_diff.min()-0.1, g_diff.max()+0.2))
p.outline_line_color = None
p.yaxis.minor_tick_line_color = None
p.xaxis.minor_tick_line_color = None
p.grid.grid_line_color = None
slope_comp = 0.75
intercept1 = plotting_utils.first_intercept(slope_comp, x.max(), g_diff.min()) -0.5
line_scale = (y.max() - y.min())/5
around_line=0.2
for i in line_scale*np.array(range(30))+intercept1:
try:
lx, ly = plotting_utils.generate_line(intercept=i, slope=slope_comp,
bounds=(x[0]-around_line, x[1]+around_line,
g_diff.min()-around_line, g_diff.max()+around_line+0.4), point=x[1])
p.line(lx, ly, color='grey', alpha=0.3)
except:
pass
bs_slope_reps, bs_intercept_reps, _, _ = plotting_utils.draw_bs_pairs_linreg(m, g_diff, size=resample_size)
p.title.text = ('Slope: '+ str(round(slope, 2)) + ' intercept: ' + str(round(intercept, 2)) +' slope 95% CI: ' +
str([round(j, 3) for j in np.percentile(bs_slope_reps, [2.5, 97.5])]) + ' ' + str(np.sum(bs_slope_reps > 0.75)))
print("Diffusive:")
a, a_CI, b, b_CI = (round(slope, 3), [round(j, 3) for j in np.percentile(bs_slope_reps, [2.5, 97.5])], round(intercept, 3), [round(j, 3) for j in np.percentile(bs_intercept_reps, [2.5, 97.5])])
print('Slope: '+ str(a) +' slope 95% CI: ' + str(a_CI) + ' Intercept: ' + str(b) +' intercept 95% CI: ' + str(b_CI))
print(a, a_CI, b, b_CI)
x_boot = np.linspace(m.min(), m.max(), 200)
y_boot = np.outer(bs_slope_reps, x_boot) + np.stack([bs_intercept_reps]*200, axis=1)
low, high = np.percentile(y_boot, [2.5, 97.5], axis=0)
p1 = np.append(x_boot, x_boot[::-1])
p2 = np.append(low, high[::-1])
p.patch(p1, p2, alpha=0.5, color='lightgrey')
p.circle(m, g_diff, color='black', size=cs)
p.line(x, y, color='black', line_width=lw, line_cap='round')
p.output_backend='svg'
plots.append(p)
#bokeh.io.show(p)
m = df_summed['log mass (g)']
g_adv = np.log10(df_summed['g_adv']*2) #double the number for g_diff since we only measured one side of the animal, they have two of each spiracle
slope, intercept = np.polyfit(m, g_adv, deg=1)
x = np.array([m.min(), m.max()])
y = slope * x + intercept
p = bokeh.plotting.figure(plot_height=300, plot_width=400, x_range=(x[0]-0.1, x[1]+0.1), y_range=(g_adv.min()-0.2, g_adv.max()+0.5))
p.outline_line_color = None
p.yaxis.minor_tick_line_color = None
p.xaxis.minor_tick_line_color = None
p.grid.grid_line_color = None
slope_comp = 0.75
intercept1 = plotting_utils.first_intercept(slope_comp, x.max(), g_adv.min())
line_scale = (y.max() - y.min())/7
around_line=0.4
for i in line_scale*np.array(range(30))+intercept1:
try:
lx, ly = plotting_utils.generate_line(intercept=i, slope=slope_comp,
bounds=(x[0]-around_line, x[1]+around_line,
g_adv.min()-around_line, g_adv.max()+around_line+0.5), point=x[1])
p.line(lx, ly, color='grey', alpha=0.3)
except:
pass
bs_slope_reps, bs_intercept_reps, _, _ = plotting_utils.draw_bs_pairs_linreg(m, g_adv, size=resample_size)
p.title.text = ('Slope: '+ str(round(slope, 2)) + ' intercept: ' + str(round(intercept, 2)) +' slope 95% CI: ' +
str([round(j, 3) for j in np.percentile(bs_slope_reps, [2.5, 97.5])]) + ' ' + str(np.sum(bs_slope_reps < 0.75)))
print("Advective:")
a, a_CI, b, b_CI = (round(slope, 3), [round(j, 3) for j in np.percentile(bs_slope_reps, [2.5, 97.5])], round(intercept, 3), [round(j, 3) for j in np.percentile(bs_intercept_reps, [2.5, 97.5])])
print('Slope: '+ str(a) +' slope 95% CI: ' + str(a_CI) + ' Intercept: ' + str(b) +' intercept 95% CI: ' + str(b_CI))
print(a, a_CI, b, b_CI)
x_boot = np.linspace(m.min(), m.max(), 200)
y_boot = np.outer(bs_slope_reps, x_boot) + np.stack([bs_intercept_reps]*200, axis=1)
low, high = np.percentile(y_boot, [2.5, 97.5], axis=0)
p1 = np.append(x_boot, x_boot[::-1])
p2 = np.append(low, high[::-1])
p.patch(p1, p2, alpha=0.5, color='lightgrey')
p.circle(m, g_adv, color='black', size=cs)
p.line(x, y, color='black', line_width=lw, line_cap='round')
p.output_backend='svg'
plots.append(p)
bokeh.io.show(plots[0])
bokeh.io.show(plots[1])
To calculate the ∆PO2 across the spiracles needed to supply beetle metabolic demand by diffusion, the metabolic rate for a quiescent beetle at a body temperature of 25 °C of a given mass was estimated from here with the following equation: $\text{log}_{10} (\text{metabolic rate } \mu \text{W}) = 3.20 + 0.75*\mathrm{log}_{10}(\text{mass (g)}$ and the assumption of 20.7 kJ/L for O2. For O2 at 25 °C, 24.5 mol/L was also assumed. Based on here, flight metabolic rate of small insects is in the range of 8x resting, whereas it is on the order of 32x resting metabolic rate in large insects. Most scarab beetles are endothermic during flight, so flight metabolic rates of these warm beetles could double this value to 64x with a 10 °C increase in thoracic temperature if the Q10 is 2. Even this calculation may be a conservative estimate of hovering metabolic rate, as oxygen consumption rose to 90x higher than those of quiescent, 25 °C fig beetles in one study, and maximal flight metabolic rates may be 1.25-2x higher than during hovering. Therefore, we estimated maximal aerobic metabolic rate during flight as 90x those of quiescent beetles. The required PO2 gradient across the spiracles to support gas exchange by diffusion at rest and during flight was calculated by rearranging our previous $G_{\text{diff}}$ equation and performing unit conversions as follows:
With this relation in hand, we plot the required oxygen partial pressure gradient across just the spiracle needed to supply the metabolic demand of the insect.
m = df_summed['log mass (g)']
g_diff = np.log10(df_summed['g_diff']*2)
slope, intercept = np.polyfit(m, g_diff, deg=1)
po2 = np.log10((1*((10**(3.20 + 0.75*m))*(1/20.7)*(1/24.5))/(df_summed['g_diff']*2)))
print(slope)
N = 200
y_top = 90
x = np.linspace(m.min(), m.max(), N)
y = np.linspace(1, y_top, N)
im = np.zeros((N, N))
for j, yi in zip(range(N), y):
for i, xi in zip(range(N), x):
im[i, j] = (yi*((10**(3.20 + 0.75*xi))*(1/20.7)*(1/24.5))/(10**(slope * xi + intercept)))
N = 800
y_top = 90
x = np.linspace(m.min(), m.max(), N)
y = np.linspace(1, y_top, N)
x_pos = []
for j, yi in zip(range(N), y):
[x_pos.append([xj, yj, zj]) for xj, yj, zj in zip(x, np.log10(yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/(10**(slope * x + intercept))),
yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/(10**(slope * x + intercept)))]
x_min, y_min, v_min = (np.array(x_pos)[:, 0].min(), np.array(x_pos)[:, 1].min(), np.array(x_pos)[:, 2].min())
x_stride, y_stride = ((np.array(x_pos)[:, 0].max() - x_min)/N, (np.array(x_pos)[:, 1].max() - y_min)/N)
im = np.ones((N, N))*v_min
for xj, yj, vj in x_pos:
im[int(np.ceil((yj-y_min)/y_stride-1)), int(np.ceil((xj-x_min)/x_stride-1))] = vj
im = scipy.ndimage.gaussian_filter(im, sigma=1)
p = bokeh.plotting.figure(tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")], plot_height=300, plot_width=400)
#p.x_range.range_padding = p.y_range.range_padding = 0
cmap = bokeh.models.LinearColorMapper(palette='Viridis256', low=im.min(), high=im.max())
cmap_low, cmap_high = (np.min(1*(10**po2)), np.max(90*(10**po2)))
cmap = bokeh.models.LinearColorMapper(palette='Viridis256', low=cmap_low, high=cmap_high)
p.image(image=[im], x=x.min(), y=np.log10(im.min()), dw=x.max()-x.min(), dh=np.log10(im.max())-np.log10(im.min()), color_mapper=cmap, level="image", )
color_bar = bokeh.models.ColorBar(color_mapper=cmap, location=(0,0), ticker=bokeh.models.BasicTicker(desired_num_ticks=12, base=10))
p.add_layout(color_bar, 'right')
p.grid.grid_line_color = None
p.patch([x.min(), x.max(), x.max(), x.min()], [np.log10((8*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min(),
np.log10((8*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
np.log10(im.min()), np.log10(im.min())], color='white', line_width=2)
p.patch([x.min(), x.max(), x.max(), x.min()], [np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min(),
np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
np.log10(im.max()), np.log10(im.max())], color='white', line_width=2)
p.line([x.min(), x.max()], [np.log10(21), np.log10(21)], color='lightgrey', line_width=2, alpha=0.75)
#p.patch([x.min(), x.max(), x.max(), x.min()], [np.log10((1*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min(),
# np.log10((1*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
# np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
# np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min()], color='white', line_width=0, alpha=0.2)
line_color='black'
#[p.line(x, np.log10((yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))), color=line_color) for yi in np.linspace(1, y_top, 10)]
[p.line(x, np.log10((yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))), color=line_color, line_width=lw, line_cap='round') for yi in [1, 8, 90]]
#[p.line(x, np.log10((yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(0.33 * x + intercept))), color='white', line_dash='dashed', line_width=5, line_cap='round') for yi in [1, 8, 90]]
dot_color='white'
dot_line='black'
dot_width = 1
cmapper = bokeh.transform.linear_cmap('c', palette='Viridis256', low=np.min(1*(10**po2)), high=np.max(90*(10**po2)))
for mult in [1, 8, 90]:
source = bokeh.models.ColumnDataSource(data=dict(x=m, y=np.log10(mult*(10**po2)), c=mult*(10**po2)))
p.circle('x', 'y', source = source, size=cs, line_color=dot_line, line_width = dot_width, fill_color=cmapper)
p.outline_line_color = None
p.yaxis.minor_tick_line_color = None
p.xaxis.minor_tick_line_color = None
p.output_backend='svg'
plots.append(p)
bokeh.io.show(bokeh.layouts.gridplot([plots[0], plots[2], plots[1]], ncols=3))
We can see that, at 90X resting and a metabolic scaling slope of $\mathrm{mass}^{0.75}$, that the partial pressure required to get enough air just across the spiracles exceeds that atmospheric oxygen partial pressure. Note that the color bar uses a linear color mapper, so the colors correspond to the colors seen the main figure panel, but the y axis of the plot does not correspond to y axis of the color bar.
In addition to exploring what the required partial pressure gradient of oxygen across the spiracles would be for flight metabolic rates of 8X or 90X resting, we also calculated these partial pressure values assuming aerobic scope (AS) for flight was 90X resting, and different scaling slopes for the metabolic rate while flying. In particular, we used the equation
where AS is the aerobic scope (90 for flight, 1 for resting) and EXP is the scaling exponent. For reasons described in our main manuscript, we used 0.67 as a conservative lower possible value flight metabolic rate scaling and 1.19 as an upper value. We make some plots of what regressions using these scaling values look like.
m = df_summed['log mass (g)']
g_diff = np.log10(df_summed['g_diff']*2)
slope, intercept = np.polyfit(m, g_diff, deg=1)
AS = 1
po2_75_resting = np.log10((1*((10**(np.log10(AS) + 3.20 + 0.75*m))*(1/20.7)*(1/24.5))/(df_summed['g_diff']*2)))
AS = 90
po2_67 = np.log10((1*((10**(np.log10(AS) +3.20 + 0.67*m))*(1/20.7)*(1/24.5))/(df_summed['g_diff']*2)))
po2_75 = np.log10((1*((10**(np.log10(AS) +3.20 + 0.75*m))*(1/20.7)*(1/24.5))/(df_summed['g_diff']*2)))
po2_119 = np.log10((1*((10**(np.log10(AS) +3.20 + 1.19*m))*(1/20.7)*(1/24.5))/(df_summed['g_diff']*2)))
def plot_bootstrap(x, y, p = bokeh.plotting.figure(), resample_size=10_000, point_color='black', band_color='lightgrey', line_color='black', lw=2, cs=12):
slope, intercept = np.polyfit(x, y, deg=1)
x_line = np.array([x.min(), x.max()])
y_line = slope * x_line + intercept
bs_slope_reps, bs_intercept_reps, _, _ = plotting_utils.draw_bs_pairs_linreg(x, y, size=resample_size)
a, a_CI, b, b_CI = (round(slope, 3), [round(j, 3) for j in np.percentile(bs_slope_reps, [2.5, 97.5])], round(intercept, 3), [round(j, 3) for j in np.percentile(bs_intercept_reps, [2.5, 97.5])])
print('Slope: '+ str(a) +' slope 95% CI: ' + str(a_CI) + ' Intercept: ' + str(b) +' intercept 95% CI: ' + str(b_CI))
print(a, a_CI, b, b_CI)
x_boot = np.linspace(m.min(), m.max(), 200)
y_boot = np.outer(bs_slope_reps, x_boot) + np.stack([bs_intercept_reps]*200, axis=1)
low, high = np.percentile(y_boot, [2.5, 97.5], axis=0)
p1 = np.append(x_boot, x_boot[::-1])
p2 = np.append(low, high[::-1])
p.patch(p1, p2, alpha=0.5, color=band_color)
p.circle(x, y, color=point_color, size=cs)
p.line(x_line, y_line, color=line_color, line_width=lw, line_cap='round')
p.output_backend='svg'
return(p)
p = bokeh.plotting.figure(plot_height=300, plot_width=400)
p.line([np.min(m), np.max(m)], [np.log10(21), np.log10(21)], color='grey', line_width=4, line_cap='round')
p = plot_bootstrap(m, po2_75_resting, point_color=bokeh.palettes.viridis(3)[0], band_color='lightgrey', line_color=bokeh.palettes.viridis(3)[0], p=p)
p = plot_bootstrap(m, po2_67, p = p, point_color=bokeh.palettes.viridis(3)[1], band_color='lightgrey', line_color=bokeh.palettes.viridis(3)[1])
#p = plot_bootstrap(m, po2_75, p=p, point_color=bokeh.palettes.viridis(3)[1], band_color='lightgrey', line_color=bokeh.palettes.viridis(3)[1])
p = plot_bootstrap(m, po2_119, p=p, point_color=bokeh.palettes.viridis(3)[2], band_color='lightgrey', line_color=bokeh.palettes.viridis(3)[2])
p.outline_line_color = None
p.yaxis.minor_tick_line_color = None
p.xaxis.minor_tick_line_color = None
p.grid.grid_line_color = None
p.output_backend='svg'
bokeh.io.show(p)
With either value for the slope, we see that large beetles would need a partial pressure difference of greater than atmospheric oxygen levels (21 kPa shown as the grey line) to get enough oxygen just across their spiracle (not to mention the rest of the tracheal system) to serve the metabolic demand.
We will now make a table with the the required partial pressure gradient with a 90x resting metabolic demand (with the various scaling exponents mentioned above) as well as some of the relevant calculated respiratory parameters.
po2_df = pd.DataFrame({'species':df_summed['species'],
'beetle mass (g)':10**m,
'summed area (mm^2)':df_summed['area']*2,
'average depth (mm)':df_summed['average depth'],
'summed G_diff (nmol/(sec*kPa))':df_summed['g_diff']*2,
'summed G_adv (cm^3/(sec*kPa))':df_summed['g_adv']*2,
'PO2_resting_0.75_exp (kPa)':10**po2_75_resting,
'PO2_flight_0.67_exp (kPa)': 10**po2_67,
'PO2_flight_1.19_exp (kPa)': 10**po2_119 })
po2_df
#po2_df.to_clipboard()
These values for $\Delta PO_2$ for a flying animal exceed the partial pressure of air in the atmosphere, both for the conservative metabolic rate scaling exponent of 0.67 or the more likely realistic exponent of 1.19! No way a large beetle can use diffusion for gas exchange, it must use active advective processes.
#TO EXPORT THE PLOTS FROM ABOVE
#bokeh.io.export_svgs(plots[2], filename="./ptest.svg")
%reload_ext watermark
%watermark -p bokeh