Branching principles of animal and plant networks identified by combining extensive data, machine learning, and modeling

TL;DR

This study combines extensive data, machine learning, and modeling to classify and understand the principles of branching networks in animals and plants, revealing key ratios that distinguish different biological systems.

Contribution

It introduces a novel classification method for diverse branching networks by integrating machine learning with new theoretical insights relating vascular form to metabolic function.

Findings

Ratios of limb radii effectively distinguish different branching networks.

Vascular and branching geometry variation persists despite convergent metabolic relationships.

A new framework links vascular form to metabolic function across species.

Abstract

Branching in vascular networks and in overall organismic form is one of the most common and ancient features of multicellular plants, fungi, and animals. By combining machine-learning techniques with new theory that relates vascular form to metabolic function, we enable novel classification of diverse branching networks--mouse lung, human head and torso, angiosperm and gymnosperm plants. We find that ratios of limb radii--which dictate essential biologic functions related to resource transport and supply--are best at distinguishing branching networks. We also show how variation in vascular and branching geometry persists despite observing a convergent relationship across organisms for how metabolic rate depends on body mass.

Figures

Tables

Table 1. Table 1 : | | | Basic measures of the different vascular branching networks studied.

Vascular branching network	Trunk radius (cm)	Trunk length (cm)	Mean tip radius (mm)	Mean tip length (mm)	Number of tips	Number of generations	Number of junctions	Total number of vessels
ML (N=1)	0.0686	0.103	0.098 (0.055)	0.709 (0.434)	688	9	660	1348
HHT (N=18)	0.383 (0.09)	2.77 (1.74)	0.855 (0.474)	6.950 (7.18)	50 (30)	6 (1)	48 (30)	1891
Balsa (N=1)	18.8	1170	5.97 (4.99)	125 (151)	357	8	292	649
Piñon (N=1)	5.73	5.4	2.30 (1.61)	24.5 (18.2)	1286	10	813	2099
Ponderosa (N=5)	2.06 (0.769)	30.6 (30.1)	2.25 (0.714)	77.8 (62.0)	31 (21)	5 (1)	23 (21)	312
Roots (N=314)	0.307 (0.293)	25.9 (23.8)	1.19 (0.838)	89.2 (74.1)	2 (2)	1 (1)	1 (1)	1231
AS/GS Tips (N=31)	0.320 (0.103)	10.8 (6.3)	0.516 (0.318)	57.8 (53.7)	15 (11)	4 (1)	12 (9)	914

Table 2. Table 2 : | | | Global scores and effect sizes for different machine learning methods and feature spaces in classifying mammal and plant datasets.

	$(r,l)$	${(r,l)}^{†}$	$({\sigma }_1,{\sigma }_2)$	($\beta ,\gamma$)	$(\overline{\beta },\overline{\gamma })$	$(\Delta \beta ,\Delta \gamma )$	$(\overline{\beta },\Delta \beta )$	$(\overline{\gamma },\Delta \gamma )$	$(\overline{\beta },\overline{\gamma },\Delta \beta ,\Delta \gamma )$
LogReg	0.82	0.54	0.52	0.52	0.59	0.55	0.59	0.53	0.58
SVM	0.88	0.56	0.52	0.57	0.62	0.59	0.64	0.58	0.67
KDE	2 $\times {10}^5$****	0.049**	0	0.065	0.11	0.31	0.72***	0.13*	0.68****

Table 3. Table 1 : | | | Results of pairwise classification of plants and animal groups using asymmetric scale factors related to hydraulics and fluid transport ( β ¯ , Δ ​ β ¯ 𝛽 Δ 𝛽 \bar{\beta},\Delta\beta ) or the strictly symmetric branching scale factors ( β , γ 𝛽 𝛾 \beta,\gamma )

	HHT*	Mouse Lung	Balsa	Piñon	Ponderosa*	GS Tips*	AS Tips*	Roots
HHT*	$-$	$2\times {10}^{-1}$	$6\times {10}^{-13}$	$1\times {10}^{-4}$	$3\times {10}^{-2}$	$9\times {10}^{-9}$	$1\times {10}^{-1}$	$2\times {10}^{-2}$
Mouse Lung	$6\times {10}^{-2}$	$-$	$5\times {10}^{-18}$	$7\times {10}^{-3}$	$7\times {10}^{-5}$	$2\times {10}^{-19}$	$2\times {10}^{-2}$	$2\times {10}^{-6}$
Balsa	$1\times {10}^{-1}$	$4\times {10}^{-4}$	$-$	$1\times {10}^{-16}$	$6\times {10}^{-2}$	$1\times {10}^{-6}$	$6\times {10}^{-2}$	$6\times {10}^{-13}$
Piñon	$3\times {10}^{-7}$	$3\times {10}^{-6}$	$2\times {10}^{-7}$	$-$	$4\times {10}^{-13}$	$3\times {10}^{-56}$	$2\times {10}^{-5}$	$2\times {10}^{-13}$
Ponderosa*	$2\times {10}^{-1}$	$7\times {10}^{-2}$	$2\times {10}^{-3}$	$4\times {10}^{-4}$	$-$	$9\times {10}^{-5}$	$3\times {10}^{-1}$	$1\times {10}^{-4}$
GS Tips*	$9\times {10}^{-3}$	$3\times {10}^{-2}$	$8\times {10}^{-4}$	$1\times {10}^{-3}$	$1\times {10}^{-1}$	$-$	$7\times {10}^{-12}$	$2\times {10}^{-11}$
AS Tips*	$2\times {10}^{-1}$	$1\times {10}^{-1}$	$7\times {10}^{-3}$	$2\times {10}^{-2}$	$4\times {10}^{-1}$	$3\times {10}^{-2}$	$-$	$5\times {10}^{-3}$
Roots	$1\times {10}^{-1}$	$6\times {10}^{-3}$	$2\times {10}^{-3}$	$2\times {10}^{-5}$	$1\times {10}^{-1}$	$3\times {10}^{-2}$	$3\times {10}^{-1}$	$-$

Table 4. Table 2 : | | | Results of pairwise classification of plant and animal groups using asymmetric scale factors that are related to space-filling only ( γ ¯ , Δ ​ γ ¯ 𝛾 Δ 𝛾 \bar{\gamma},\Delta\gamma ) or space-filling and hydraulics ( β ¯ , Δ ​ β , γ ¯ , Δ ​ γ ¯ 𝛽 Δ 𝛽 ¯ 𝛾 Δ 𝛾 \bar{\beta},\Delta\beta,\bar{\gamma},\Delta\gamma ).

	HHT*	Mouse Lung	Balsa	Piñon	Ponderosa*	GS Tips*	AS Tips*	Roots
HHT*	$-$	$7\times {10}^{-2}$	$8\times {10}^{-1}$	$6\times {10}^{-3}$	$1\times {10}^{-1}$	$1\times {10}^{-2}$	$3\times {10}^{-2}$	$8\times {10}^{-3}$
Mouse Lung	$1\times {10}^{-1}$	$-$	$1\times {10}^{-1}$	$2\times {10}^{-3}$	$4\times {10}^{-1}$	$4\times {10}^{-1}$	$9\times {10}^{-2}$	$9\times {10}^{-6}$
Balsa	$9\times {10}^{-5}$	$2\times {10}^{-5}$	$-$	$2\times {10}^{-2}$	$1\times {10}^{-2}$	$1\times {10}^{-3}$	$1\times {10}^{-2}$	$3\times {10}^{-2}$
Piñon	$4\times {10}^{-11}$	$4\times {10}^{-8}$	$1\times {10}^{-8}$	$-$	$3\times {10}^{-1}$	$6\times {10}^{-2}$	$5\times {10}^{-1}$	$3\times {10}^{-1}$
Ponderosa*	$8\times {10}^{-3}$	$9\times {10}^{-3}$	$7\times {10}^{-3}$	$6\times {10}^{-5}$	$-$	$3\times {10}^{-1}$	$5\times {10}^{-1}$	$1\times {10}^{-1}$
GS Tips*	$2\times {10}^{-4}$	$6\times {10}^{-4}$	$1\times {10}^{-3}$	$1\times {10}^{-5}$	$2\times {10}^{-1}$	$-$	$3\times {10}^{-1}$	$9\times {10}^{-4}$
AS Tips*	$5\times {10}^{-2}$	$2\times {10}^{-2}$	$1\times {10}^{-1}$	$1\times {10}^{-2}$	$4\times {10}^{-1}$	$1\times {10}^{-1}$	$-$	$3\times {10}^{-1}$
Roots	$3\times {10}^{-4}$	$1\times {10}^{-4}$	$8\times {10}^{-8}$	$1\times {10}^{-3}$	$2\times {10}^{-2}$	$1\times {10}^{-3}$	$4\times {10}^{-2}$	$-$

Table 5. Table 3 : | | | Results of pairwise classification of plant and animal groups using asymmetric scale factors that are related to hydraulics and space filling and that are associated with the central means of asymmetric branching ( β ¯ , γ ¯ ¯ 𝛽 ¯ 𝛾 \bar{\beta},\bar{\gamma} ) or scale factors associated with variation in asymmetric branching ( Δ ​ β , Δ ​ γ Δ 𝛽 Δ 𝛾 \Delta\beta,\Delta\gamma ).

	HHT*	Mouse Lung	Balsa	Piñon	Ponderosa*	GS Tips*	AS Tips*	Roots
HHT*	$-$	$4\times {10}^{-1}$	$4\times {10}^{-3}$	$5\times {10}^{-3}$	$3\times {10}^{-1}$	$3\times {10}^{-1}$	$2\times {10}^{-1}$	$4\times {10}^{-2}$
Mouse Lung	$2\times {10}^{-1}$	$-$	$3\times {10}^{-2}$	$1\times {10}^{-3}$	$3\times {10}^{-1}$	$3\times {10}^{-1}$	$3\times {10}^{-1}$	$7\times {10}^{-3}$
Balsa	$3\times {10}^{-2}$	$8\times {10}^{-7}$	$-$	$2\times {10}^{-4}$	$2\times {10}^{-1}$	$8\times {10}^{-2}$	$7\times {10}^{-2}$	$1\times {10}^{-5}$
Piñon	$4\times {10}^{-2}$	$1\times {10}^{-2}$	$3\times {10}^{-6}$	$-$	$1\times {10}^{-1}$	$8\times {10}^{-2}$	$2\times {10}^{-1}$	$8\times {10}^{-2}$
Ponderosa*	$3\times {10}^{-2}$	$1\times {10}^{-3}$	$3\times {10}^{-2}$	$1\times {10}^{-4}$	$-$	$6\times {10}^{-1}$	$5\times {10}^{-1}$	$9\times {10}^{-2}$
GS Tips*	$3\times {10}^{-6}$	$1\times {10}^{-11}$	$3\times {10}^{-7}$	$2\times {10}^{-12}$	$9\times {10}^{-2}$	$-$	$5\times {10}^{-1}$	$3\times {10}^{-2}$
AS Tips*	$1\times {10}^{-1}$	$1\times {10}^{-2}$	$9\times {10}^{-2}$	$6\times {10}^{-2}$	$3\times {10}^{-1}$	$4\times {10}^{-3}$	$-$	$1\times {10}^{-1}$
Roots	$1\times {10}^{-1}$	$4\times {10}^{-3}$	$2\times {10}^{-2}$	$1\times {10}^{-1}$	$9\times {10}^{-2}$	$1\times {10}^{-3}$	$4\times {10}^{-1}$	$-$

Table 6. Table 4 : | | | Results of pairwise classification of Ponderosa individuals using scale factors that are related to hydraulics, and fluid transport ( β ¯ , Δ ​ β ¯ 𝛽 Δ 𝛽 \bar{\beta},\Delta\beta ) and that are associated with strictly symmetric branching ( β , γ 𝛽 𝛾 \beta,\gamma ) or asymmetric branching.

	POND03	POND05	POND06	POND07	POND16
POND03	$-$	$5\times {10}^{-2}$	$5\times {10}^{-1}$	$7\times {10}^{-1}$	$3\times {10}^{-1}$
POND05	$5\times {10}^{-1}$	$-$	$3\times {10}^{-2}$	$8\times {10}^{-2}$	$3\times {10}^{-1}$
POND06	$3\times {10}^{-1}$	$4\times {10}^{-1}$	$-$	$4\times {10}^{-1}$	$3\times {10}^{-1}$
POND07	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$5\times {10}^{-1}$	$-$	$6\times {10}^{-1}$
POND16	$2\times {10}^{-1}$	$2\times {10}^{-1}$	$2\times {10}^{-1}$	$2\times {10}^{-1}$	$-$

Table 7. Table 5 : | | | Results of pairwise classification of the human head and torso individuals using scale factors that are related to hydraulics and fluid transport and that are associated with strictly symmetric branching ( β , γ 𝛽 𝛾 \beta,\gamma ) or asymmetric branching ( β ¯ , Δ ​ β ¯ 𝛽 Δ 𝛽 \bar{\beta},\Delta\beta ).

	HHT01	HHT02	HHT03	HHT04	HHT05	HHT06	HHT07	HHT08	HHT09	HHT10	HHT11	HHT12	HHT13	HHT14	HHT15	HHT16	HHT17	HHT18
HHT01	$-$	$6\times {10}^{-1}$	$7\times {10}^{-1}$	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-2}$	$5\times {10}^{-1}$	$6\times {10}^{-1}$	$5\times {10}^{-1}$	$1\times {10}^{-1}$	$5\times {10}^{-1}$	$7\times {10}^{-1}$	$3\times {10}^{-1}$	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$6\times {10}^{-1}$	$6\times {10}^{-1}$
HHT02	$3\times {10}^{-1}$	$-$	$7\times {10}^{-1}$	$7\times {10}^{-1}$	$7\times {10}^{-1}$	$3\times {10}^{-1}$	$6\times {10}^{-2}$	$7\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-2}$	$5\times {10}^{-1}$	$6\times {10}^{-1}$	$2\times {10}^{-1}$	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$7\times {10}^{-1}$	$7\times {10}^{-1}$
HHT03	$5\times {10}^{-1}$	$3\times {10}^{-1}$	$-$	$8\times {10}^{-1}$	$6\times {10}^{-1}$	$3\times {10}^{-1}$	$5\times {10}^{-2}$	$7\times {10}^{-1}$	$7\times {10}^{-1}$	$7\times {10}^{-1}$	$1\times {10}^{-1}$	$6\times {10}^{-1}$	$7\times {10}^{-1}$	$2\times {10}^{-1}$	$4\times {10}^{-1}$	$7\times {10}^{-1}$	$7\times {10}^{-1}$	$7\times {10}^{-1}$
HHT04	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$4\times {10}^{-1}$	$-$	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-2}$	$7\times {10}^{-1}$	$6\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-2}$	$6\times {10}^{-1}$	$5\times {10}^{-1}$	$2\times {10}^{-1}$	$3\times {10}^{-1}$	$6\times {10}^{-1}$	$7\times {10}^{-1}$	$7\times {10}^{-1}$
HHT05	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$3\times {10}^{-1}$	$4\times {10}^{-1}$	$-$	$2\times {10}^{-1}$	$9\times {10}^{-3}$	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$4\times {10}^{-1}$	$7\times {10}^{-3}$	$2\times {10}^{-1}$	$4\times {10}^{-1}$	$1\times {10}^{-1}$	$3\times {10}^{-1}$	$6\times {10}^{-1}$	$6\times {10}^{-1}$	$7\times {10}^{-1}$
HHT06	$7\times {10}^{-2}$	$3\times {10}^{-2}$	$1\times {10}^{-1}$	$6\times {10}^{-2}$	$2\times {10}^{-2}$	$-$	$5\times {10}^{-1}$	$3\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-1}$	$6\times {10}^{-1}$	$9\times {10}^{-1}$	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$2\times {10}^{-1}$	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$3\times {10}^{-1}$
HHT07	$3\times {10}^{-1}$	$2\times {10}^{-2}$	$2\times {10}^{-1}$	$2\times {10}^{-1}$	$1\times {10}^{-2}$	$2\times {10}^{-1}$	$-$	$1\times {10}^{-2}$	$2\times {10}^{-1}$	$3\times {10}^{-1}$	$6\times {10}^{-1}$	$2\times {10}^{-1}$	$3\times {10}^{-3}$	$2\times {10}^{-1}$	$3\times {10}^{-1}$	$1\times {10}^{-2}$	$2\times {10}^{-1}$	$7\times {10}^{-2}$
HHT08	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$4\times {10}^{-1}$	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-2}$	$2\times {10}^{-1}$	$-$	$3\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-2}$	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-1}$	$1\times {10}^{-1}$	$5\times {10}^{-1}$	$8\times {10}^{-1}$	$8\times {10}^{-1}$
HHT09	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$5\times {10}^{-1}$	$6\times {10}^{-1}$	$3\times {10}^{-1}$	$2\times {10}^{-1}$	$4\times {10}^{-1}$	$5\times {10}^{-1}$	$-$	$3\times {10}^{-1}$	$4\times {10}^{-1}$	$7\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$6\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$6\times {10}^{-1}$
HHT10	$4\times {10}^{-1}$	$3\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-1}$	$2\times {10}^{-1}$	$2\times {10}^{-1}$	$2\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-1}$	$-$	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$5\times {10}^{-1}$	$1\times {10}^{-1}$	$3\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$
HHT11	$3\times {10}^{-1}$	$3\times {10}^{-2}$	$2\times {10}^{-1}$	$1\times {10}^{-1}$	$1\times {10}^{-2}$	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$2\times {10}^{-1}$	$5\times {10}^{-1}$	$2\times {10}^{-1}$	$-$	$6\times {10}^{-1}$	$3\times {10}^{-2}$	$7\times {10}^{-1}$	$4\times {10}^{-1}$	$5\times {10}^{-2}$	$2\times {10}^{-1}$	$5\times {10}^{-2}$
HHT12	$5\times {10}^{-1}$	$3\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$2\times {10}^{-1}$	$1\times {10}^{-1}$	$2\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$6\times {10}^{-1}$	$3\times {10}^{-1}$	$-$	$3\times {10}^{-1}$	$8\times {10}^{-1}$	$3\times {10}^{-1}$	$4\times {10}^{-1}$	$8\times {10}^{-1}$	$6\times {10}^{-1}$
HHT13	$6\times {10}^{-1}$	$5\times {10}^{-1}$	$6\times {10}^{-1}$	$7\times {10}^{-1}$	$4\times {10}^{-1}$	$4\times {10}^{-2}$	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$6\times {10}^{-1}$	$3\times {10}^{-1}$	$1\times {10}^{-1}$	$4\times {10}^{-1}$	$-$	$1\times {10}^{-1}$	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$5\times {10}^{-1}$
HHT14	$2\times {10}^{-1}$	$2\times {10}^{-2}$	$2\times {10}^{-1}$	$1\times {10}^{-1}$	$3\times {10}^{-2}$	$1\times {10}^{-1}$	$2\times {10}^{-1}$	$2\times {10}^{-1}$	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$3\times {10}^{-1}$	$5\times {10}^{-1}$	$1\times {10}^{-1}$	$-$	$4\times {10}^{-1}$	$2\times {10}^{-1}$	$5\times {10}^{-1}$	$3\times {10}^{-1}$
HHT15	$5\times {10}^{-1}$	$2\times {10}^{-1}$	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$2\times {10}^{-1}$	$6\times {10}^{-2}$	$3\times {10}^{-1}$	$5\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-1}$	$4\times {10}^{-1}$	$5\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-1}$	$-$	$2\times {10}^{-1}$	$3\times {10}^{-1}$	$3\times {10}^{-1}$
HHT16	$5\times {10}^{-1}$	$6\times {10}^{-1}$	$4\times {10}^{-1}$	$6\times {10}^{-1}$	$5\times {10}^{-1}$	$2\times {10}^{-2}$	$5\times {10}^{-2}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$3\times {10}^{-1}$	$6\times {10}^{-2}$	$3\times {10}^{-1}$	$6\times {10}^{-1}$	$6\times {10}^{-2}$	$4\times {10}^{-1}$	$-$	$5\times {10}^{-1}$	$7\times {10}^{-1}$
HHT17	$4\times {10}^{-1}$	$6\times {10}^{-1}$	$6\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$1\times {10}^{-1}$	$3\times {10}^{-1}$	$6\times {10}^{-1}$	$5\times {10}^{-1}$	$3\times {10}^{-1}$	$4\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$2\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$-$	$7\times {10}^{-1}$
HHT18	$5\times {10}^{-1}$	$3\times {10}^{-1}$	$5\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-1}$	$9\times {10}^{-2}$	$1\times {10}^{-1}$	$5\times {10}^{-1}$	$3\times {10}^{-1}$	$7\times {10}^{-1}$	$1\times {10}^{-1}$	$5\times {10}^{-1}$	$4\times {10}^{-1}$	$4\times {10}^{-1}$	$3\times {10}^{-1}$	$4\times {10}^{-1}$	$4\times {10}^{-1}$	$-$

Table 8. Table 6 : | | | Results of pairwise classification of tree tips species using scale factors associated with strictly symmetric branching ( β , γ 𝛽 𝛾 \beta,\gamma ), or asymmetric branching, hydraulics, and fluid transport ( β ¯ , Δ ​ β ¯ 𝛽 Δ 𝛽 \bar{\beta},\Delta\beta ).

	Maple	Oak	Robinia	Whitefir	Dougfir	Whitepine
Maple	$-$	$3\times {10}^{-2}$	$7\times {10}^{-1}$	$4\times {10}^{-10}$	$1\times {10}^{-2}$	$8\times {10}^{-3}$
Oak	$4\times {10}^{-1}$	$-$	$3\times {10}^{-2}$	$9\times {10}^{-6}$	$4\times {10}^{-3}$	$1\times {10}^{-1}$
Robinia	$5\times {10}^{-1}$	$2\times {10}^{-1}$	$-$	$3\times {10}^{-6}$	$4\times {10}^{-2}$	$3\times {10}^{-2}$
Whitefir	$3\times {10}^{-1}$	$1\times {10}^{-1}$	$4\times {10}^{-1}$	$-$	$2\times {10}^{-1}$	$3\times {10}^{-1}$
Dougfir	$4\times {10}^{-1}$	$2\times {10}^{-1}$	$5\times {10}^{-1}$	$5\times {10}^{-1}$	$-$	$1\times {10}^{-1}$
Whitepine	$5\times {10}^{-1}$	$3\times {10}^{-1}$	$6\times {10}^{-1}$	$6\times {10}^{-1}$	$6\times {10}^{-1}$	$-$

Equations

$$\bar{\beta }=\frac{{\beta }_1+{\beta }_2}{2}\bar{\gamma }=\frac{{\gamma }_1+{\gamma }_2}{2}$$

$$\Delta \beta =\frac{{\beta }_1-{\beta }_2}{2}\Delta \gamma =\frac{{\gamma }_1-{\gamma }_2}{2}$$

$$(\bar{\beta }+\Delta \beta {)}^2+(\bar{\beta }-\Delta \beta {)}^2=1$$

$$(\bar{\gamma }+\Delta \gamma {)}^3+(\bar{\gamma }-\Delta \gamma {)}^3=1$$

$${\sigma }_1=\frac{\ln (r_{c1}/r_p)}{\ln (l_{c1}/l_p)}{\sigma }_2=\frac{\ln (r_{c2}/r_p)}{\ln (l_{c2}/l_p)}$$

$${\sigma }_1=\frac{\ln (\bar{\beta }+\Delta \beta )}{\ln (\bar{\gamma }+\Delta \gamma )}{\sigma }_2=\frac{\ln (\bar{\beta }-\Delta \beta )}{\ln (\bar{\gamma }-\Delta \gamma )}$$

$$\theta =\frac{\ln (2^N)}{\ln (2^N)+\ln (1-{\nu }^{N+1})-\ln ({\nu }^N(1-\nu ))}$$

$$\theta \approx \frac{\ln (2)}{\ln (2)-\ln (2{\bar{\beta }}^2\bar{\gamma })-\ln (1+\frac{2\Delta \beta \Delta \gamma }{\bar{\beta }\bar{\gamma }}+\frac{\Delta {\beta }^2}{{\bar{\beta }}^2})}$$

$$TPR$$

$$FPR$$

$$B=B_0{\left(\frac{M}{M_0}\right)}^{\theta }$$

$$B_{\text{textsc}TIP}{N}_{\text{textsc}TIPS}=B_0{\left(\frac{M}{M_0}\right)}^{\theta }$$

$$N_{\text{textsc}TIPS}={\left(\frac{V_{\text{textsc}TOT}}{V_0}\right)}^{\theta }$$

$$V_{\text{textsc}TOT}=V_{\text{textsc}N,\text{textsc}TOT}\sum _{j=0}^N{\nu }^{-j}$$

$$V_{\text{textsc}TOT}=V_{\text{textsc}N,\text{textsc}TOT}\frac{1-{\nu }^{-(N+1)}}{1-{\nu }^{-1}}$$

$$V_{\text{textsc}TOT}=\{\begin{array}{ll}{V}_{\text{textsc}N,\text{textsc}TOT}\frac{1-{\nu }^{-(N+1)}}{1-{\nu }^{-1}}& \text{mbox}for\nu \ne 1\\ V_{\text{textsc}N,\text{textsc}TOT}\frac{N+1}{{\nu }^N}& \text{mbox}for\nu \approx 1\end{array}$$

$$\theta =\{\begin{array}{ll}\frac{\ln (2^N)}{\ln (2^N)+\ln (1-{\nu }^{N+1})-\ln ({\nu }^N(1-\nu ))}& \text{mbox}for\nu \ne 1\\ \frac{\ln (2^N)}{\ln ((N+1)2^N)-N\ln (\nu )}& \text{mbox}for\nu \approx 1\end{array}$$

$$\ln (B)=\ln (B_0)+\theta \ln \left(\frac{M}{M_0}\right)$$

$$\frac{\partial }{\partial \ln (M)}=\frac{1}{\ln (2)}\frac{\partial }{\partial N}$$

$$\frac{{\partial }^2\ln (B)}{\partial (\ln (M){)}^2}=\frac{2}{\ln (2)}\frac{\partial \theta }{\partial N}+\frac{N}{\ln (2)}\frac{{\partial }^2\theta }{\partial N^2}$$

$$\frac{{\partial }^2\ln (B)}{\partial (\ln (M){)}^2}=\frac{1}{x^4}\left\{2x^3+2N^{2}x{\left(\frac{\partial x}{\partial N}\right)}^2-4N{x}^2\frac{\partial x}{\partial N}-N^2{x}^2\frac{{\partial }^{2}x}{\partial N^2}\right\}$$

Full text

Title: Branching principles of animal and plant networks identified by combining extensive data, machine learning, and modeling

Short Title: Functional groups of biologic branching

Authors: Alexander Byers Brummer1,2,3∗,

Panagiotis Lymperopoulos1, Jocelyn Shen4,

Elif Tekin1,2, Lisa P. Bentley5, Vanessa Buzzard6,

Andrew Gray6, Imma Oliveras7,

Brian J. Enquist6,8, Van M. Savage1,2,3,8

1Institute for Quantitative and Computational Biology, University of California, Los Angeles

2Department of Computational Medicine, David Geffen School of MedicineUniversity of California, Los Angeles

3Department of Ecology and Evolutionary Biology, University of California, Los Angeles

4 Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology

5 Department of Biology, Sonoma State University

6 Department of Ecology and Evolutionary Biology, University of Arizona

7 Environmental Change Institute, University of Oxford

8 Santa Fe Institute

∗To whom correspondence should be addressed; E-mail: abrummer@ucla.edu

One Sentence Summary: Vascular theory and machine-learning enable principled classification of the form and function of biologic branching in mammals and plants.

Abstract: Branching in vascular networks and in overall organismic form is one of the most common and ancient features of multicellular plants, fungi, and animals. By combining machine-learning techniques with new theory that relates vascular form to metabolic function, we enable novel classification of diverse branching networks—mouse lung, human head and torso, angiosperm and gymnosperm plants. We find that ratios of limb radii—which dictate essential biologic functions related to resource transport and supply—are best at distinguishing branching networks. We also show how variation in vascular and branching geometry persists despite observing a convergent relationship across organisms for how metabolic rate depends on body mass.

Introduction

It is a great challenge to decipher which features of biological branching networks are shared, which are different, and when these differences matter [1, 2]. For instance, branching in plant and animal networks exhibit strikingly similar features despite profound physiological and environmental differences (e.g., carbon dioxide and sap versus oxygen and blood, mobile versus stationary organisms, heart and pulsatile flow versus non-pulsatile flow) [3, 4, 5, 6, 7, 8, 9, 10, 11]. Similarly, differences in loopiness and “noisiness” are well documented between vascular branching in tumors or stroke-damaged tissue versus healthy tissue [2, 12, 13]. The shared branching features are argued to lead to functional convergence in plant and animal networks via biological rates despite the notable physiological differences just listed [14, 15, 16]. Yet, the extent of shared versus distinct branching features has not been systematically and quantitatively analyzed across plants and animals in the same study. Consequently, there is a need to understand the forces that shape the full spectrum of form and function in branching networks (Figure 1a).

The classification of branching architectures is historically based on coarse qualitative differences in morphological features. Examples include: classifying lobes of the liver based on independent blood supply [17]; or the paired/un-paired ordering of plant leaves along a stem [18]. Recent efforts have identified gene expression profiles related to branching phenotypes, with examples of branching in the developing lung as being planar versus tetrahedral in orientation [19], or branching in the developing kidney based on the number of terminal vessels downstream from (or distal to) two sibling branches [20]. However, these empirically motivated classifications still fall short of relating patterns in vascular form to biophysical and biomechanical function.

With recent advances in automated methods of image analysis developed by us and others [4, 5, 21], increasing amounts of data are becoming available to tackle these problems. The tools that are missing are efficient and accurate algorithms for categorizing branching across whole networks and different organisms. In this paper we apply machine-learning methods to theoretically-informed feature spaces to leverage all available information and technology to achieve these goals.

We analyze the largest-ever compilation of branching network data, with over 58 distinct networks and approximately 8,000 vessels or tree limbs (see Table 1). We collected these data over the last decade for both mammalian cardiovascular systems and plant architecture in both angiosperms and gymnosperms. Two mammalian networks are studied, the first being the major arterial branching junctions of the human head and torso (HHT) for 18 adult individuals (H. sapiens) collected using contrast-enhanced magnetic resonance angiography on a 3 Tesla Siemens Trio scanner with voxel dimensions between $700\times 700\times 800\mu m^{3}$ and $800\times 800\times 900\mu m^{3}$ [5]. The second mammalian network is the full pulmonary vascular branching of one wild-type adult mouse lung (ML) (M. musculus) collected using a combination of vascular casting with MICROFIL and micro computed tomography on a $\mu$CT 40, ScanCo Medical scanner with 10$\mu$m isotropic voxel spacing [9]. All mammalian network data was acquired using the open source software Angicart [22].

The plant networks consist of: 1. whole, above-ground, adult trees for one Balsa (O. pyramidale), one Piñon (P. edulis), and five Ponderosa pines (P. ponderosa) [3], 2. an array of angiosperm root clusters belonging to Andean tropical montane cloud forests [23], and 3. a collection of 50 cm long clippings of the terminal ends of canopy branches from three species each of angiosperms (AS Tips) and gymnosperms (GS Tips) comprised of Maple (A. grandidentatum), Scrub Oak (Q. gambelii), Robinia (R. neomexicana), White Fir (A. concolor), Douglas Fir (P. menziesii), and White Pine (P. strobiformis). Tree measurements—all done destructively by hand—are of the external branching structures (limbs), not the xylem that are directly responsible for water transport. Scaling relationships for the external limbs directly determine similar relationships for the internal xylem based on previous empirical studies [24, 25] and established branching theory [11, 26], thus enabling comparisons of plant and animal networks for the structure, flow, and function in the present study [3, 21, 27].

To search for patterns, machine-learning is often applied to the full set of untransformed, standardized raw data. These raw data thus represent one feature space, yet there are always infinitely more choices of feature spaces based on specific combinations, subsets, mathematical operations (e.g., logarithms or ratios), or other transformations of the raw data (Figure 1c, d, e). Informed choices of feature space hold the promise of greatly improving the convergence time, accuracy, and inference of machine-learning algorithms. Here we show how crucial this choice can be and the roles that our understanding of the underlying biology can play in its selection. We further demonstrate that this approach identifies key strengths and weaknesses in the theory used to guide the transformations, and thus informs our understanding, or lack thereof, of the underlying biology and physics.

The default choice for feature spaces for our networks would be the centered and standardized raw data—all vessel radii and lengths for branching networks. However, theory grounded in evolution, biology, and physics predicts that the parent-to-child ratios of radii and length—along with associated scaling exponents throughout the networks [10, 11]—encapsulate the most biologically-informative properties because they are directly tied to organismic function. Specifically, numerous models tie these ratios to the ability of branching networks to efficiently fill space and to deliver resources [7, 10, 11, 26, 28]. The fine-scale relationships between fluid flow, global vascular or branching architecture, and vessel or branch morphology are indeed complex [29]. Despite this, much information can be gleaned from the connections between the radial scale factors and hydrodynamics and the length scale factors and space-filling as first-order effects [1, 3, 5, 9]. As candidates for second-order effects, we also examine branch slenderness exponents. These couple the radial and length scale factors and inform the likelihood that a branch will experience gravitational buckling under its own weight [6, 27, 30].

We use recent theory developed by some of us (Brummer, et al. in [28]) for the asymmetric branching patterns that are pervasive throughout our data. In this theory the two sibling vessels (labelled $c1$ and $c2$, Figure 1b) and the parent vessel (labelled $p$) are combined to give two radial scale factors $\beta_{1}=r_{c1}/r_{p}$ and $\beta_{2}=r_{c2}/r_{p}$ and two length scale factors $\gamma_{1}=l_{c1}/l_{p}$ and $\gamma_{2}=l_{c2}/l_{p}$. Thus the average radial and length scale factors are

[TABLE]

To capture sibling branch asymmetry the difference radial and length scale factors are

[TABLE]

Corresponding constraint equations for area-preserving and space-filling branching—used in canonical optimization models—are

[TABLE]

Separately, susceptibility to gravitational buckling is quantified in the slenderness exponents [30], which relate the scaling of radii to the scaling of lengths as

[TABLE]

It is not a priori obvious which combinations of the scale factors will work best as a feature space for discriminating vascular networks. If dynamics of blood flow dominate the formation and evolution of vascular architecture, then variation in scale factors involving vessel radius would be expected to be most informative because vascular theory and empirical evidence show blood flow is most strongly determined by vessel radius. For example, it is well documented that as blood flow transitions from pulsatile to non-pulsatile, so too does the scaling of vessel radii from the squared scaling (scaling exponent $=2$) of Eq. (3) to cubic scaling (scaling exponent $=3$) similar to Eq. (4) [10]. This quantitative shift would then show up as a difference between classified groups in our data that should be detectable by, and informative to, our machine learning algorithms. If the space-filling constraints and body plan of the organism primarily determine vascular architecture, then variation in scale factors for vessel lengths should best discriminate. Moreover, average properties might be shared across species while differences or variation around these average properties could reflect distinct selective pressures that can be used to discriminate types of networks and branching principles. Alternatively, some selective pressures could change the average properties yet share the same values of variation and asymmetry. Finally, if resilience to gravitational buckling determines branching form then the slenderness exponents should differentiate between those organisms susceptible to buckling (plants) and those that are not (mammals).

To test and quantify all of these possibilities, we generate distributions of our data for combinations of the raw and standardized radius and length measurements $(r,l)$ and $(r,l)^{\dagger}$ (where ${\dagger}$ represents the centered and standardized radii and lengths), the slenderness exponents ($\sigma_{1},\sigma_{2}$), and of the symmetric and asymmetric scale factors ($\beta_{1}$, $\beta_{2}$, $\gamma_{1}$, $\gamma_{2}$, $\bar{\beta}$, $\bar{\gamma}$, $\Delta\beta$, $\Delta\gamma$) for the combined mammal and plant networks. We first examine the performance of several standard machine-learning techniques to categorize our network data [31, 32]. We use principle components analysis (PCA) to examine feature space variance (Figure 1f-g), and compare the results of support vector machine (SVM), logistic regression (LR), and kernel density estimation (KDE) machine-learning methods (Table 2 and Figure 2a-c). Uncertainty is controlled for by graphing the rates of true positive detection versus false positive detection in a one-versus-all comparison between the different classifiers being used while varying the significance of classification (Figure 2). See Supplementary Materials for additional detail on training and testing protocol. Upon identifying which method has the greatest overall classification success, we then examine which regions in the plant and mammal feature space drive classification and correspond to different species or tissues (Figures 3-4). Here we account for uncertainty by bootstrapping-with-replacement on the training and testing groups when examining better method at a fixed level of classification significance. Finally, by drawing on metabolic scaling theory—the prediction that the scaling of organism metabolism with mass is determined by vascular geometry—we examine how these different feature spaces constrain variation in estimates of the scaling exponent for organismal metabolic rate (Figure 5).

Results

We demonstrate the importance of choosing theoretically-informed feature spaces over raw data to classify vascular organisms relating form to function. Classification using only raw data (branch radii and lengths) results only in size-based categorization, an approach that can distinguish between a mouse lung and a Balsa tree, but is not easily applicable to similarly sized organisms or tissues (Figure 1c). Once networks are normalized for size, distributions of the raw data are greatly overlapped [3, 5] (Figure 1d) and machine-learning methods applied to the raw data cannot distinguish the networks (Table 2). We thus conclude that our theoretically-informed feature spaces are objectively superior at categorizing branching networks over raw data. In addition, these theoretically-informed feature spaces facilitate much easier translation into known biological principles and constraints on biologic function related to blood flow, space-filling, and metabolic rate.

Importantly, not all theory-motivated features improve classification. Table 2 shows that the slenderness scaling exponents $(\sigma_{1},\sigma_{2})$ do no better than random chance. Two likely explanations are that the definitions in Eq. (5) simultaneously de-couple branching asymmetry and couple radial and length scaling. This is supported by Figure 1f-g, where PCA loadings between $\Delta\beta,\sigma_{1}$, and $\sigma_{2}$ are either directly correlated (PCA 2), or are linear combinations of each other (PCA 3 and 4). Table 2 shows that radial scaling with asymmetry ($\bar{\beta},\Delta\beta$) outperforms length scaling with asymmetry ($\bar{\gamma},\Delta\gamma$). Thus, it appears that transformations that suppress asymmetry and enhance length—such as the slenderness scaling exponents—act to obscure the defining features between the mammal and plant data considered.

Differences in machine learning performance across methods (KDE, SVM, LR) is due to the multivariate structure of the feature spaces being studied and the chosen machine learning method. In particular, the KDE method excels at resolving the multimodality [33] that characterizes the radial scale factors for the plant dataset (Figures 2a and 3a-b). Since the distribution means are approximately equivalent, the SVM and LR methods are strongly influenced by outliers and the higher moments comprising the mammal dataset (Figures 2b) [32].

Comparing across all methods and features spaces, we find that the combination of the KDE method and the average and difference scale factors for radius ($\bar{\beta},\Delta\beta$) are the most effective for classifying branching network data (Figure 2, Table 1, and Supplementary Material Tables S1-S3). The fact that variation in the feature space for radial scale factors is the best has strong implications about what functional features form the major distinctions between biologic networks. Specifically, our empirical finding of the primacy of information based on scaling ratios of radii strongly suggests that hydrodynamic principles are the primary drivers of vascular branching patterns and overall network form.

Multiple theories of vascular networks, as well as basic physics and fluid mechanics, dictate that rates of fluid flow are largely governed by the total cross-sectional area of vessels or limbs, which can be exactly related to the ratios of scaling radii used in our feature space [7]. Importantly, theory recently developed by us demonstrates that there can exist a range of morphologies that still adhere to these area-preserving—pulsatile flow in mammals or external branching in plants—or area-increasing predictions—non-pulsatile flow in mammals [28]. In contrast, variation in the ratios of vessel lengths is more strongly tied to the ability of the vascular network to fill the body. Thus, length ratios appear to either encapsulate little important information about the differences among biologic networks, or they may not adequately capture the key properties of space filling for the architecture of vascular networks [3, 5, 28, 34]. Having identified the best performing feature space and machine learning method, we now delve deeper into the variation in the architecture and functional properties of vascular networks.

Focusing on the KDE method we see that mammalian branching exhibits more area-increasing branching than plants (Figure 3b(i)). Area-increasing branching is necessary to simultaneously increase total surface area for oxygen and metabolite transport and to slow blood flow as it travels from the heart to the capillaries and transitions from pulsatile to non-pulsatile flow, the latter phenomenon not present in plants. However, values of $\bar{\beta}\approx 1.0$ and $\Delta\beta\approx 0$ represent a deviation from the theoretical predictions of $\Delta\beta=0$ and $\bar{\beta}=1/2^{1/3}\approx 0.794$ for the non-pulsatile flow expected in this region. This marked increase in cross-sectional area is shared by both the HHT and ML networks as indicated by the nearly null relative abundances of these two networks (Figure 3c(i)) as well as by the insignificant $p$-value score of 0.2 from the global-level implementation of the KDE method (See Supplementary Table S1). This suggests that transitions in blood flow type from pulsatile to non-pulsatile may occur across a greater range of branching generations, and begin nearer to the heart, than in current theory [7, 16].

The majority of plant networks adhere to area-preservation while exhibiting a greater tendency than mammals to branch asymmetrically (specifically the Balsa, Piñon, Ponderosas, and GS Tips, Figure 3a and b). Within the plants we find that differentiation is driven at the species level (Figure 3c(ii)-(iv) and Supplementary Table S1), unrelated to plant categorization as angiosperm or gymnosperm. For example the Balsa, an angiosperm, is the only species present in both the positive and negative asymmetry types (Figures 1b) as demonstrated by being the only network with its standard deviation outside the null expectation in Figure 3c(ii) and (iv). Thus, the Balsa consists of two unique branching motifs that distinguish it from GS Tips, and the Piñon and Roots that have large relative abundances in one region each—the negative asymmetric branching of motif c(iv) and the symmetric branching of motif c(iii), respectively.

Mechanisms for the asymmetry and motifs observed in the plant radial scale factors are likely due to functional trait plasticity associated with light-seeking behavior, self- and wind-induced pruning, gap-filling, and other environmental stressors. However, making quantitative connections remains an open challenge [6, 27, 30]. For example, the slenderness scaling exponents can be calculated for all six scenarios in Figure 3c(ii)-(iv), beginning with the first generation of child branches. Expressing Eq. (5) in terms of the average and difference scale factors

[TABLE]

Thus, the slenderness scaling exponents for each plant motif are: $\sigma_{1}=-12.3,\sigma_{2}=1.4$ for motif (ii); $\sigma_{1}=-19.0,\sigma_{2}=5.1$ for motif (iii); and $\sigma_{1}=-4.8,\sigma_{2}=0.54$ for motif (iv). Biomechanical theory that applies columnar (Euler beam) buckling to branching systems demonstrates that slenderness exponents of $\sigma\geq 1$ are structurally advantageous for plant architectures as they push the locations of breakage points into the canopy as opposed to the trunk. Yet, the slenderness exponents we calculate for the observed motifs do not entirely agree with this framework, despite adherence of the radial scale factors in motifs (ii) and (iv) to the area-preserving branching constraint of Eq. (3). This disparity may lie in two sources: (i) the slenderness exponent formula of Eq. (5) was originally derived using symmetrically branched networks and (ii) the length scaling exponents involved in Eq. (5). The latter issue we now investigate.

Connecting length-based categorization to mechanism—the space-filling constraint of Eq. (4) —remains a challenge. The combination of the KDE method and length scale factor feature space ($\bar{\gamma},\Delta\gamma$) identified only one region of significance. In this region differentiation is driven by the plants, specifically the Piñon and Roots (Figure 4 and Supplementary Table S3). This is despite the large amount of variance explained by the length scale factors in the principle components analysis (Figure 1). The single region driving differentiation corresponds to average length scale factors $\bar{\gamma}<1$. This effect would normally result in an increase in the slenderness exponent, Eq. (5), driving gravitationally induced buckling (self-pruning) to occur in the canopy instead of at the trunk ($\sigma\geq 1$) [30]. However, median values of $\sigma$ for both the plants and mammals were approximately 0.2, far below the needed theoretical threshold of $\sigma=1$. We interpret this deviation from expected biomechanics as an indicator that the length scale factors, as defined, are poor features for characterizing vascular or branching architecture.

The inability of the length scale factors to inform classification between networks suggests several scenarios. Two contrasting and extreme scenarios are that either a universal architecture or a completely random architecture is being followed by both the mammals and plants [9]. This result is unlike the radial scaling that is strongly coupled to hydraulics. Current theory suggests that the architecture associated with length scaling is guided by the principles of space-filling fractals [7, 10, 11, 28]. However, large deviations are observed between the joint distributions of the length scale factors and the theoretical curves determined by the space-filling conservation equation (Figure 4a). A third scenario is that there exists a disconnect between how length scale factors are conventionally defined in simplified models versus how they are measured in complicated natural systems. All three scenarios support the need for including missing constraints, variables, and assumptions (e.g. branching angles, multi-fractal scaling, etc.), or alternative mathematical frameworks [34, 36, 37, 38, 39].

To better understand the physiological and biological implications of these categorizations, we examine the influence of asymmetric branching on estimates of biological rates—specifically the metabolic scaling exponent $\theta$ that canonically relates metabolic rate $B$ to body mass $M$ as $B\propto M^{\theta}$. Previous studies spanning orders of magnitude in body mass have shown that $\theta$ converges on a value near 3/4, yet exhibits variation specific to mammals or plants [14, 15, 16].

To probe this variation we use branching data to estimate metabolic scaling (Figure 5) by directly accounting for network geometry and size [7, 10, 28],

[TABLE]

where $N$ is the total number of branching generations in the network and $\nu$ represents volumetric scaling—the ratio of the sum of the volumes of both child branches to the volume of the parent branch. Specification of $\nu$ allows estimation of $\theta$ under different model assumptions for symmetric ($\nu=2\beta^{2}\gamma$) or asymmetric ($\nu=2\bar{\beta}^{2}\bar{\gamma}+4\bar{\beta}\Delta\beta\Delta\gamma+2\bar{\gamma}\Delta\beta^{2}$) branching. We also use a regression method between the number of terminal branches $N_{{\textsc{\tiny TIPS}}}$ and total volume $V_{{\textsc{\tiny TOT}}}$ distal to a given branch ($N_{{\textsc{\tiny TIPS}}}\propto V_{{\textsc{\tiny TOT}}}^{\theta}$) that does not depend directly on geometry (See Supplementary Materials).

We find that asymmetric branching increases the predicted values of metabolic scaling exponents when compared to the symmetric- and regression-based methods (Figure 5a). This is due to all networks exhibiting some length asymmetry, and more importantly suggests that previous studies have underestimated metabolic scaling exponents by not accounting for such variation [3, 40, 41].

To understand which different scale factors are primarily responsible for observed variation in the predicted metabolic scaling exponents we focus on the asymmetric version of Eq. (7). Estimated metabolic scaling exponents are graphed for each individual organism in terms of the average scale factors ($\bar{\beta},\bar{\gamma}$) in Figure 5b and difference scale factors ($\Delta\beta,\Delta\gamma$) in Figure 5c. We compare these graphs against the corresponding theoretical predictions reproduced from Brummer et al. in [28] where we have graphed the approximate form of Eq. (7),

[TABLE]

assuming small volume scaling ($\nu<1$), generationally large networks ($N>>1$), and enforcing area-preserving and space-filling (Eqs. (3) and (4)).

We observe a striking amount of grouping among the mammals and plants when graphing the metabolic scaling exponent $\theta$ versus the average radial and length scale factors $\bar{\beta}$ and $\bar{\gamma}$ (Figure 5b). This indicates that, of all the features and data considered, the average scale factors ($\bar{\beta}$, and $\bar{\gamma}$) are the primary determinants of variation in the metabolic scaling exponent and thus organism function.

In contrast to previous theory and importantly for understanding how diverse branching architectures could lead to universal scaling exponents, we find near constancy of the metabolic scaling exponent despite large fluctuations in length scaling (Figure 5c). These shared exponents are likely driven by the little to no radial asymmetry observed in mammalian networks and suggests that variation in length asymmetry ($\Delta\gamma$) in vascular networks has little influence on whole organism metabolic function in the presence of symmetric radial branching ($\Delta\beta=0$).

Figures 5b-c demonstrate marked deviation in the observed grouping (or lack thereof) between the empirically based predictions of metabolic scaling from Eq. (7) and the constraint-based theoretical predictions of metabolic scaling from Eq. (8). To explore this deviation we calculate curvature between metabolic rate and mass in log-log space (Supplementary Materials). When branching networks are strictly assumed to be very large, ($N>>1$) and decreasing in volume in all segments across any generation ($\nu<1$, Eq. (8)), we predict zero curvature, regardless of the extent of branching asymmetry. When accounting for variation in network size and volume scaling (Eq. (7)) we predict positive (concave up) curvature (Figure 5d). These predictions are both in agreement with respiration-based studies of mammals [16], and demonstrate the need for theories of metabolic scaling that incorporate the finite-size of the network. Furthermore, we predict that curvature decreases to zero with increasing network size, or generation $N$, in agreement with respiration-based studies of plants [15]. These results can be informative for future studies that simultaneously connect branching patterns and vascular data to ontogenetic- and size-based shifts in organismal metabolism. Such shifts are observed in growth and reproduction curves for tumors [42], plants [15], mammals [43], and fish [44].

Discussion

Machine learning is a powerful tool, but often considered a black box. We show that by combining machine learning with mechanistic theory it can be made more effective and provide insight into physiological mechanism. Here we take a first step towards building that bridge by using mechanistic theory of vascular networks to choose better feature spaces. In so doing, we achieve a two-fold, mutually reinforcing benefit: 1. We achieve better results for categorizing networks than if we used either the raw feature space or a mechanistically inspired feature space that predicts only one morphology (the symmetric scale factors $\beta$ and $\gamma$), and 2. Results are much more interpretable. For example, the best-performing feature spaces—the asymmetric ratios of vessel and limb radii, $\bar{\beta}$ and $\Delta\beta$—are explicitly and naturally tied to specific mechanisms—hydrodynamic constraints and resource flow—and allow for variation in form while still following these constraints. Additionally, the under-performing feature spaces—ratios of vessel and limb lengths, $\bar{\gamma}$ and $\Delta\gamma$—identify what may be potential weaknesses in the theory and avenues for new development to provide greater specificity. Alternatively, the inability of the ratios of vessel and limb lengths to classify between mammals and plants may be pointing to broadly shared architectural principles that are not specific to mammals or plants. This is despite different mechanistic demands, such as structural support for plants.

The results of this study also serve to inform our understanding of the physiological pressures that determine convergence in organismal form and function. We find that variation in vascular based estimates of metabolic scaling exponents—in particular curvature—is primarily determined by variation in the average scale factors ($\bar{\beta}$ and $\bar{\gamma}$), symmetric radial branching, and relative network size. This result helps to resolve some of the contradictory size-based observations in variation in metabolic scaling between mammals and plants [15, 16]. It emphasizes the needs to develop models of vascular networks that can better account for finite-size measurements (clustered sampling versus whole network measurement), to acquire comprehensive datasets that span the entirety of the vascular branching structures being studied, and to simultaneously acquire respiration based measurements of organismal metabolism.

In this direction, a shortcoming of our model of metabolic curvature is its complete inability to capture negative curvature. This scenario arises when examining individual growth curves in mammals, plants, and tumors when the metabolic scaling exponent decreases from linear to sub-linear (typically from 1 to 3/4) [15, 42, 43]. Failure to capture this essential biologic feature should spur continued development of theories of metabolic scaling and vascular branching.

Finally, this work has implications for several fields, spanning bio-mechanical and physiological imaging and theory to machine learning and biomedical applications. Incorporating topological features—connectivity and loops—and branching angles could enhance categorization methods because these features provide structural integrity and redundancy to damage in plant leaves and in capillaries [45, 46, 47]. Additional measures that capture organ and organismal physiology could provide further insight and tests. Examples for mammalian tissues include flow reserve—the change in blood flow between normal and dilated vessel states—or blood perfused for a given vascular tree [48, 49]. New applications of tomographic imaging and computer vision techniques to plants—light detection and ranging (LiDAR) and positron emission tomograpahy (PET)—are greatly expanding digitized plant architecture datasets and allowing for the direct inclusion of branch angles and xylem and phloem transport measurements as part of the biological feature space [21, 41, 50, 51]. Simultaneously, advances in medical imaging and vascular segmentation algorithms are leading to datasets of fully connected branching and blood vessel networks [52]. Such expansive datasets previously unavailable will allow for comprehensive testing of vascular branching theories where, in principle, machine learning based motif identification could be used to digitally regenerate branching networks using iterated function systems [36, 39].

Utilizing more robust applications of machine-learning methods and model complexity might help improve classification based on raw data and should improve classification using feature spaces based on theory as well. In closing, this work provides a proof-of-principle that a mechanistically-based automatic classification and detection scheme for vascular networks could have application in medical diagnostics for long term progressive disease (e.g., tumor growth). Here, classification is driven by outlier detection between vascular networks surrounding and comprising tumors compared against verified healthy vascular neteworks [53]. Such an application would serve as a new dimension in radiomic studies where the detection and classification of tumors based on vascular branching is wholly absent [2, 54] and could provide an alternative measure of tumor growth and development [55].

Funding

This work was supported by the National Science Foundation Grant 1254159

Author Contributions

ABB, BJE, and VMS conceived the analysis, ET, LPB, VB, AG, and IO acquired the data. ABB, PL, and JS, analyzed the data. ABB and VMS developed theory. ABB, BJE, and VMS wrote the manuscript. All authors edited the manuscript.

Competing Interests

The authors declare that they have no competing financial interests.

Correspondence

Correspondence and requests for materials should be addressed to ABB (email: abrummer@ucla.edu).

Data Availability

The datasets and code that support the findings of this study have been uploaded as part of the supplementary material.

Supplementary materials

Supplementary Text

Figures S1 to S3

Tables S1 to S6

References 1-19

Supplementary Datasets

Supplementary Code

SUPPLEMENTARY INFORMATION

Contents

1. A Classifier methods

2. B Feature space selection

3. B.1 Raw radius and length measurements

4. B.2 Scale factor feature spaces

5. C Normalizing results from classifier methods

6. C.1 Kernel density estimation

7. C.2 Logistic regression and support vector machine

8. C.3 ROC comparison of methods

9. C.4 Data grouping

10. D Derivation of exact metabolic scaling exponent formula

11. E Derivation of curvature in metabolic rate versus mass

12. F Kernel Density Estimation Results

Appendix A Classifier methods

To generate probability models using the support vector machine and logistic regression methods, approximately 75% of each combined plant and and combined mammal dataset was taken as the training sample and the remaining 25% used as the testing sample. This was done after randomization within respective datasets to minimize the chances of accidentally removing an entire individual or species. Then, using the radial scale factor feature space ($\bar{\beta}$ and $\Delta\beta$) the two groups of mammals and plants were compared using three methods: a support vector machine (SVM) model and a logistic regression (LR) model, and a non-parametric kernel density estimator (KDE) method.

The support vector machine method, a supervised machine learning algorithm for classification problems, plots data points in an n-dimensional space and draws a decision boundary by maximizing the margin between points from different classes. To generate a support vector machine model for our datasets, we used the SVC (support vector classification) function from the Python scikit-learn package. When running SVC, we used the polynomial kernel with a degree of 2.

The LR method generates probabilities of how likely certain data points are going to fall under a certain class, based on the logit function $1/(1+e^{-x})$. To produce a non-linear LR model for our two-dimensional datasets, we used the Python scikit-learn logistic regression function and added three nonlinearity columns, $x_{1}^{2},x_{2}^{2}$, and $x_{1}\times x_{2}$. All pairwise combinations were run on this LR model and the probabilities were recorded.

For both the SVM method and the LR method, a polynomial kernel was used. Through a series of tests between different kernels, including linear and radial basis function, the polynomial kernel yielded the highest accuracy score when classifying the testing data. The LR and SVM methods differ in their nonlinearities, where the LR method uses an added nonlinearity term and the SVM method uses a radial basis function kernel. The probabilities for the LR method are assigned using the logit function, and the probabilities assigned using the SVM method are based on the point distances from the decision boundary. For both of these methods, testing points are classified based on the value of their score. When using training data that is equally split between the two categories then points that receive scores of 0.5 or greater are classified as one group, and scores below 0.5 are classified as the other group.

The kernel density estimator (KDE) technique introduced by Duong et al. [1, 2] generates non-parametric, multi-dimensional probability distributions, $\mathbb{P}_{i}(\textbf{x})$, of vascular traits, x, for each testing group, $i=A,B$. These distributions are then compared against one another for their extent of uniqueness, or non-overlap, represented by the test statistic $T=\int[\mathbb{P}_{A}(\textbf{x})-\mathbb{P}_{B}(\textbf{x})]^{2}d\textbf{x}$. This test can be thought of as a multi-dimensional generalization of the two-sample Kolmogorov-Smirnov (KS) test, where significance of classification is conventionally communicated through p-values. We use the nominal threshold of $p\leq 0.05$ as a threshold for significance when using the KDE method.

The KDE method can be applied both globally [1] and locally [2]. At the global level, the test statistic $T$ is calculated as described above (and in more detail in [1]), and converted to a $p$-value using standard tables. While $p=0.05$ is used as the nominal threshold for significance, it should be noted that $p$ values ranging orders of magnitude in size from $p=10\times 10^{-2}$ to $p=10\times 10^{-20}$ are observed, and thus interpreted as relative levels of significance. The local application of the KDE method is effectively an inversion of the calculation for the test statistic $T$. Here, one sets the desired threshold for significance, or the value of the test statistic $T$, then calculates which regions in the vascular trait-space, x correspond to the chosen test statistic.

Appendix B Feature space selection

Here we describe our methods for selecting the feature space. We tested each machine learning method on a variety of vascular trait feature spaces. Specifically, we tested the feature spaces of: raw and standardized diameter and length measurements; the slenderness scaling exponents, $\sigma_{1}$, and $\sigma_{2}$; the symmetric WBE diameter and length scale factors, $\beta$ and $\gamma$; and five combinations of the asymmetric scale factors. The five combinations of asymmetric scale factors were: the average scale factors $\bar{\beta}$ and $\bar{\gamma}$; the difference scale factors $\Delta\beta$ and $\Delta\gamma$; the radial scale factors $\bar{\beta}$ and $\Delta\beta$; the length scale factors $\bar{\gamma}$ and $\Delta\gamma$; and all four asymmetric scale factors $\bar{\beta},\bar{\gamma},\Delta\beta,$ and $\Delta\gamma$. These results are presented in Table 2. We also conducted a principle components analysis (PCA) on the asymmetric scale factors to identify which combinations of scale factors explain the most variance in the datasets (Figure 1f-g). Performing a PCA on all vascular variables is non-trivial due to the non-random presence of empty-cells in the dataset [3].

B.1 Raw radius and length measurements

Prior to transforming the radial and length measurements of our vascular datasets as motivated by scaling theory [4, 5, 6], we applied all three of the logistic regression (LR), support vector machine (SVM), and kernel density estimator (KDE) methods on the raw, untransformed data. This was done in two ways, first on the data after all measurements were converted to meters (Figure 1c), then again after performing a standardized transformation by translating each species- or group-level distribution to be centered about zero, and then normalizing by the respective standard deviations (Figure 1d).

Classifying vascular data based on metric size is both trivial and uninformative. Small networks (mouse lung) are clearly distinct from large networks (whole trees), and varying degrees of overlap will exist in the intermediate range of all other networks (Figure 1c). All three methods yield significant global classification scores, as demonstrated by the example scores of LR: 82%, SVM: 88%, and KDE: T-stat $\approx 2\times 10^{5}$ found in Table 2. One can certainly argue that classifying networks in this manner is possible, even though it is simply demonstrating that networks of differing size are distinguishable. However, this approach does not provide an obvious path toward understanding at a mechanistic level why certain patterns are observed beyond the simple size-based classification. Furthermore, when applying these methods for the purposes of disease detection one is oftentimes examining healthy and diseased tissues that are of similar size classes. Thus, the utility of size-based classification is rendered moot, and we must transform.

The most common transformation is to standardize the data such that the distribution means are centered about the origin and to normalize by the variance [7, 8] (Figure 1d). While this approach has the benefit of removing the global size-based hierarchy between the networks, it fails to address the common pattern of local size-based hierarchy that is commonly found within networks [4, 5, 9]. Specifically, the abundance-size distribution of vessels in a vascular network is approximately exponential due to the fact that the vessels, on average, decrease in size at every bifurcation. Thus, classification becomes immediately obscured, as demonstrated by the noticeably decreased example scores of LR: 54%, SVM: 56%, and KDE: T-stat $=0.049$ (with a corresponding $p$-value less than 0.001) in Table 2. The fact that the KDE method still yielded a significant score is a consequence of this method excelling at detecting multimodality between distributions. Yet, even though the KDE method does yield a significant score, the effect size is negligible and we still have the original problem of how to interpret the results, only now it is further compounded after having performed the standardization transformation. It is possible to consider alternative transformations based on mechanical principles—flow-rates and pressures [10]—yet the problem of the hierarchy of sizes returns, only now with different physical units. As a consequence of these complications in analyzing raw data, we turn to classification based on scale factor variables as guided by the metabolic scaling theory literature.

B.2 Scale factor feature spaces

To investigate the scale factor feature spaces, all three machine learning methods were tested on different combinations of candidate feature spaces, and a principle component analysis (PCA) was conducted on the asymmetric scale factor feature space. Results from the application of different machine learning methods are presented in Table 2 and graphs from the PCA are presented in Figure 1f-g. Table 2 shows that the top two combinations of features for classifying animal and plant vascular networks are the set of all asymmetric scale factors ($\bar{\beta},\Delta\beta,\bar{\gamma},\Delta\gamma$) and the radial scale factors ($\bar{\beta},\Delta\beta$) . That the set of all asymmetric scale factors performs the best is not surprising as this is the most comprehensive set of data for the vascular networks considered. However, it does not help to identify which subset of variables are primarily responsible for classification over other variables. To do that, we focus our attention to the PCA results.

A common feature of vascular network data is the large variation observed in the the scaling of branch lengths [5, 4, 9], often times exhibited over many orders of magnitude. In Figure 1f we see that the first principle component explains 34.2% of the variance in the data, and is composed primarily of the length scale factors ($\bar{\gamma},\Delta\gamma$). Interestingly, the length scale factors are not powerful for classification purposes (see Table 2 and Figure 4 in the main text), even though they are responsible for much of the variance in the data. On the other hand, we find that correlations between the radial scale factors account for 17.3% and 16.4% of the variance through the second and third principle components. This result, in combination with the fact that the radial scale factors had the best global classification scores in Table 2, led to the selection of the radial average and difference scale factors as the final choice for the feature space.

Appendix C Normalizing results from classifier methods

While the three machine learning methods tested effectively perform the same task (e.g. classification of data), the manner in which this is done varies significantly, and thus requires careful consideration when trying to compare results. This difficulty is compounded by the inclusion of a variable classification sensitivity level. Here we describe the process used to standardize classification output from the three methods of Kernel Density Estimation (KDE), Support Vector Machine (SVM), and Logistic Regression (LR). Only one variable feature space was used for the comparison between methods, and that was the radial scale factor feature space ($\bar{\beta},\Delta\beta$). This feature space was chosen as it performed best for classification across all methods (Table 2, as well as explaining up to 33% of the variance in all of the asymmetric scale factor variables (Figure 1f-g). The testing groups used were that of mammal and plant. These were chosen to increase the performance of methods reliant on dataset size for training and testing. Finally, we present receiver operating characteristic (ROC) curves that are used to perform the final comparison once standardization has been achieved.

C.1 Kernel density estimation

The non-parametric kernel density estimation procedure put forth by Duong et al. [1, 2] tests for uniqueness and overlap between two different probability distributions generated from empirical data. Probability distributions, $\mathbb{P}_{i}(\textbf{x, X})$ are defined on a discretized feature space, where x represents the discretized coordinates in the feature space, X represents actual empirically measured values, and $i=A$ or $B$ represents our known classifier. When performing the local significant difference test, one must first set the $p$-value that denotes “significance”. Conventionally this value is set at either $p=0.01$ (as done for the analysis presented in the main text) or $p=0.05$. The KDE method then identifies regions of the feature space where the squared-difference between the probabilities, $\mathbb{P}_{i}(\textbf{x, X})$, is less than the selected significance level. Finally, by examining which probability, $\mathbb{P}_{i}(\textbf{x, X})$, is greatest within each region, one can identify which classifier is driving differentiation. Thus, one can measure the performance of the KDE method at a given significance level by counting the number of correctly and incorrectly classified points with respect to each region (Figure S1).

By varying the $p$-value one can vary the relative size of the classified regions of the test to examine the efficacy of classification (Figure S1). For the KDE method, the $p$-value was varied by integer orders of magnitude from $10^{-12}$ to $10^{4}$ to ensure that the limiting scenarios of zero points classified and all points classified were contained. Then, for each sensitivity level, correctly and incorrectly classified points were tallied for further analysis and comparison in a one-versus-all framework.

C.2 Logistic regression and support vector machine

The logistic regression (LR) and support vector machine (SVM) methods represent two approaches at using supervised machine learning methods for classification [7, 11]. These methods differ from the KDE approach by using a training set of data to partition the feature space into two separable regions (separated by the decision boundary). Then, points from a testing set of data are assigned a classification score based on their positions in the feature space with respect to the decision boundary. In Figure S1 are graphs of results from the LR method and the SVM method. In these examples, the decision boundary corresponds to a probability score of 0.5, and any testing point with a score between 0 and 0.5 is classified as a plant, while a score between 0.5 and 1 is classified as a mammal.

To compare to the output of the KDE approach, an analogue of significance to the KDE approach must be defined for the LR and SVM approaches. In the context of the LR and SVM approaches this was done by defining significance as the distance a probability score is from 0.5. Regions of equal significance were determined by binning along the probability score axis. Thus, points with probability scores in the ranges of $[0,0.05]$ for mammals or $[0.95,1]$ for plants would all be characterized as equally, highly significant predictions (or in terms of the KDE method, would correspond to the next lowest possible $p$-value). The bins were then successively enlarged to reflect a decrease in test significance. So, for mammals, the bins used were $[0,0.05]$, $[0,0.1]$, $[0,0.15]$, $[0,0.2]$, $[0,0.25]$, $[0,0.3]$, $[0,0.35]$, $[0,0.4]$, $[0,0.45]$, $[0,0.5]$. For plants, the bins used were $[0.95,1]$, $[0.9,1]$, $[0.85,1]$, $[0.8,1]$, $[0.75,1]$, $[0.7,1]$, $[0.65,1]$, $[0.6,1]$, $[0.55,1]$, $[0.5,1]$. In Figure S1 are graphs of results from the LR and SVM methods showing the successive binning approach. Finally, for each binned region corresponding to varying levels of significance, correctly and incorrectly classified points were identified for comparison in a one-versus-all framework.

C.3 ROC comparison of methods

To finally compare the three methods of kernel density estimation (KDE), logistic regression (LR), and support vector machine (SVM), we graphed receiver operating characteristic (ROC) curves of the results of the methods as sensitivities were varied [7, 11]. ROC curves are graphs of a methods true positive rate (TPR) versus its false positive rate (FPR). Due to the manner in which we are comparing machine learning methods, we technically have three classes for any given significance level: mammals, plants, and undetected. Thus, we must use a one-versus-all framework for calculating true and false positives. In calculating the TPR and FPR, we first choose which category we are “detecting”, say plants, then the TPR and FPR can be calculated as,

[TABLE]

where $N_{plant\pm}$ corresponds to the number of data points correctly (or incorrectly) classified as plant as denoted by the $+$ sign ($-$ sign) within the green plant contours in Figure 2, and $N_{not~{}plant\pm}$ corresponds to the number of data points correctly (or incorrectly) classified as not plant external to the green plant contours. A TPR and an FPR is then calculated for each level of significance (defined by a $p$-value for the KDE approach, or a probability bin for the LR or SVM approaches). Finally, the ROC curve can be graphed, as presented in Figure 2c.

Conventionally, when comparing classification schemes using an ROC curve, the best method is identified as whichever method sits most in the upper-left corner of the graph, or which ever has the greatest “area under the curve” (AUC), as this represents a maximal true positive rate and a minimal false positive rate. To use the ROC curves for the different methods of KDE, LR, and SVM, we can use Eq. (LABEL:eq:plant_positive) and its mammal classifying analogue, resulting in Figure 2c. Here we can observe that overall the KDE method outperforms the LR and SVM methods as having either a lower FPR for a given TPR (plant as positive class graph on left in Figure 2c) or a greater TPR for a given FPR (mammal as positive class graph on right in Figure 2c).

C.4 Data grouping

Once selection of the classifier was made, the subgroups of datasets were prepared for higher resolution classification. Using the KDE method, individual HHT and ponderosas were found indistinguishable, as were species within the GS Tips and AS Tips groups, hence their final merging into larger datasets (see Tables S4-6). We obtained 8 different major species/groups: HHT, mouse lung, Balsa, Piñon, Ponderosa, GS Tips, AS Tips, and Roots, each with 6 recorded variables: $\bar{\beta},\bar{\gamma},\Delta\beta,\Delta\gamma,$ and the merging of $\beta_{1}$ and $\beta_{2}$ into a single distribution as well as for $\gamma_{1}$ and $\gamma_{2}$. Specific groupings of the scale factor variables used were: ($\beta_{1},\beta_{2};\gamma_{1},\gamma_{2}$) as a two-dimensional distribution representative of the symmetric formalism; ($\bar{\beta},\bar{\gamma},\Delta\beta,\Delta\gamma$) as a four-dimensional distribution representative of the asymmetric formalism; ($\bar{\beta},\bar{\gamma}$) as a two-dimensional distribution for average scaling; ($\Delta\beta,\Delta\gamma$) as a two-dimensional distribution for difference scaling; ($\bar{\beta},\Delta\beta$) as a two-dimensional distribution for radial scaling; and ($\bar{\gamma},\Delta\gamma$) as a two-dimensional distribution for length scaling. For the full list of comparison results using the KDE method, see Tables S1-3.

Appendix D Derivation of exact metabolic scaling exponent formula

Here we present a derivation of the metabolic scaling exponent under the general assumptions of the West, Brown, Enquist model for vascular branching [12, 13, 6]. This approach deviates from conventional derivations in that zero approximations will be made regarding network size or the particulars of volume scaling outside of the bounds of self-similarity. Some discussion will be presented at the end regarding how network geometry can be included in terms of symmetric or asymmetric branching.

The typical starting point for modeling metabolic scaling is Kleiber’s Law [kleiber_hilgardia_1932], the empirically motivated, power-law relationship between organismal basal metabolic rate and mass, expressed as,

[TABLE]

Here $B$ represents the measured metabolic rate, $M$ the mass, $\theta$ the metabolic scaling exponent observed to cluster around 3/4 [14, 15], and $B_{0}$ and $M_{0}$ normalization constants. Treating metabolic rate as the combined sum of metabolism of every terminal branch in an organism, we can substitute $B=B_{\textsc{\tiny TIP}}N_{\textsc{\tiny TIPS}}$, where $B_{\textsc{\tiny TIP}}$ represents total metabolism per terminal branch and $N_{\textsc{\tiny TIPS}}$ the total number of metabolizing terminal branches. Doing so results in,

[TABLE]

Next, mass is substituted with total volume. The validity of this substitution is the result of having optimized the geometrical scaling of a hierarchically branching vascular network against the dual demands of hydrodynamic resistance to fluid flow and fractal space-filling [12, 13, 6]. Performing the substitution gives us,

[TABLE]

where we have cancelled out $B_{\textsc{\tiny TIP}}$ with $B_{0}$. At this point it is important to pause and recognize Eq. S.4 as a method by which one can estimate the metabolic scaling exponent in a vascular organism free from explicitly imposing assumptions regarding network geometry via symmetric or asymmetric branching. One powerful aspect of Eq. S.4 is that it can be applied recursively on any individual vascular branching network. Specifically, for every branch (or vessel) in a network, the total number of distal terminal branches can be represented by $N_{\textsc{\tiny TIPS}}$, and the total volume of all distal branches are represented by $V_{\textsc{\tiny TOT}}$. Thus, a standard major axis regression analysis can be performed to determine the value of $\theta$ that corresponds to an individual organism as per this model. An example of such an analysis being performed on the mouse lung data set is presented in Figure S2.

It is at this point that we must specify in greater detail the functional form of the total volume of the vascular network, $V_{\textsc{\tiny TOT}}$. Assuming that we are working with a strictly self-similar, hierarchically branching, pipe-like model, then the total volume of the network can be expressed as

[TABLE]

where $j$ and $N$ represent the $j^{th}$ and $N^{th}$ generations of the network, $V_{{\textsc{\tiny N}},{\textsc{\tiny TOT}}}$ is the total volume of the terminal ($N^{th}$) generation, and $\nu$ represents the ratio of total volume from sibling branches to their respective parent branch. For example, in a bifurcating symmetric network, $\nu=2(\pi r_{j+1}^{2}l_{j+1})/(\pi r_{j}^{2}l_{j})=2\beta^{2}\gamma$. Recognizing Eq. S.5 as a geometric series, we can write it’s exact form as,

[TABLE]

which is valid for all values of $\nu$ except for $\nu\approx 1$. It is not uncommon to find individual organisms with $\nu\approx 1$. This scenario is problematic if using the above formula due to its asymptotic nature. However, this problem can be remedied by using L’Hôpital’s rule, producing the piecewise result,

[TABLE]

Upon substituting Eq. S.7 into Eq. S.4, we arrive at the following exact, piecewise function for the metabolic scaling exponent,

[TABLE]

where we set $V_{{\textsc{\tiny N}},{\textsc{\tiny TOT}}}=N_{{\textsc{\tiny TIPS}}}V_{{\textsc{\tiny TIP}}}$, $V_{0}$ was cancelled out with $V_{{\textsc{\tiny TIP}}}$, and we have restricted ourselves to considering strictly bifurcating networks such that $N_{{\textsc{\tiny TIPS}}}=2^{N}$. One benefit of Eq. S.8 is that we can easily examine the functional relationship between metabolic rate and the scaling of volume in a general sense, as presented in Figure S3. Examining the influence of symmetric or asymmetric branching can then be done separately.

Several important features are present in Figure S3. Most notably is the effect that network size, $N$, has on varying the sensitivity of the metabolic scaling exponent to the value of the volume scale factor, $\nu$. For all values of network size, the metabolic scaling exponent is bounded between zero and one. For the case of large networks ($N=100$ in Figure S3), the metabolic scaling exponent nearly reaches the asymptotic value of $\theta=1$ when $\nu=1$, where as for smaller size networks (lesser values of $N$), the rate of increase in $\theta$ is markedly less. This result lends support to the argument that the approximate form of Eq. S.8, namely $\theta\approx\ln(2)/(\ln(2)-\ln(2\beta^{2}\gamma))$ for a symmetrically branching network, holds in the limit that $N>>1$ and $\nu<1$. This approximate form of $\theta$ is presented in Figure S3 by the solid, blue line.

A second important result present in Figure S3 is the fact that all network sizes can take on the empirically observed value of $\theta=3/4$ as long as $\nu$ is large enough. This result supports previous arguments regarding transitions in blood flow related to shifts in the scaling of radii. In these circumstances, the radial scaling transitions from cross-sectional area-preserving ($\beta^{2}=1/2$ for a symmetrically bifurcating network) to area-increasing ($\beta^{3}=1/2$). Due to the conservation fluid flow-rates across branching junctions, this latter scaling acts to slow the flow of blood. However, a simultaneous result of changing the scaling of radii in such a way is to increase the volume scale factor, assuming that there is no change in the scaling of lengths. Thus, this predicted trade-off between decreasing network size and increasing volume scaling would appear consistent with the biological demands associated with the cardiovascular system in mammals. As for plants, ontogenetic shifts in metabolic scaling (as measured through total respiration) have been reported between sapling and mature trees and shrubs, yet not in the context of the measurement of vascular branching traits [16]. Thus, it would be interesting to see if future studies corroborate trade-offs between metabolic scaling exponents, vascular branching traits, and network size.

In using Eq. S.8, one must specify (and thereby either measure or estimate) the total number of branching generations in a vascular network. Unfortunately this is not a particularly straightforward task. With the conventional WBE generational labeling scheme, the generation index is advanced once a branching occurs for all branches within the previous generation. For highly asymmetric networks this can result in having branches of wildly differing size be members of the same generation. Alternative labelling schemes have been proposed outside the context of the WBE framework, specifically that of Horton-Strahler [17, 18]. While capable of handling moderate levels of asymmetry, one can show that the Horton-Strahler labelling scheme is insufficient at handling networks resembling the fishbone branching pattern. For the purposes of this work, we adopt an approach based on the asymptotically symmetric expectation that for a network with $N$ generations, there should be approximately $2^{N}$ terminal branches. Thus, from a count of the number of terminal branches, $N_{{\textsc{\tiny TIPS}}}$, one can estimate the number of generations as $N=\ln(N_{{\textsc{\tiny TIPS}}})/\ln(2)$.

Lastly, in the main texts standard errors are reported for calculated values of the metabolic scaling exponent in Figure 4a. These errors were determined from standard methods, namely $\sigma_{\theta}^{2}\approx\left\|\frac{\partial\theta}{\partial\nu}\right\|^{2}\sigma_{\nu}$, where $\sigma_{\nu}$ is determined from the chosen symmetric or asymmetric scale factor variables.

Appendix E Derivation of curvature in metabolic rate versus mass

Here we present a derivation of the curvature in the metabolic rate, $B$, of an organism in terms versus its mass, in log-log space. Beginning again with Kleiber’s Law after having linearized the equation, we have,

[TABLE]

As we are interested in examining the mass related curvature in metabolic rate, we must evaluate $\partial^{2}\ln(B)/\partial(\ln(M))^{2}$. Since $\theta$ is effectively a mass-dependent quantity through its dependence on network size via total generation number $N$, we need to find an appropriate substitute for mass. Assuming that the total mass of a vascular organism can be expressed asymptotically as the sum total of all cells serviced by the vascular architecture, then $M\approx M_{{\textsc{\tiny TIP}}}2^{N}$. Thus, we can build a derivative operator as,

[TABLE]

Using this operator, curvature in metabolic rate can be expressed as follows,

[TABLE]

where we have maintained the explicit notation for mass, $M$, on the left hand side as a reminder of what the quantity on the right represents. Now the problem is reduced to calculating derivatives of the metabolic scaling exponent $\theta$ with respect to total number of generations $N$. An important and immediate consequence of Eq. S.11 is that the approximate ($\nu<1$) and asymptotic ($N>>1$) version of Eq. (S.8), $\theta\approx\ln(2)/\left[\ln(2)-\ln(\nu)\right]$, results in zero curvature due to its invariance in relative network size. This prediction is consistent with alternative parameterizations of curvature in metabolic scaling [19], and could inform observed deviations from previous theoretical predictions of metabolic scaling exponents as discussed in the main text.

Focusing on the version of Eq. (S.8) for $\nu\neq 1$, we find that the curvature in $B$ can be expressed finally in the relatively compact form of,

[TABLE]

where $x=\ln\left(V_{{\textsc{\tiny TOT}}}/V_{{\textsc{\tiny TIP}}}\right)$. Graphs of Eq. S.12 are presented in the main text, where we find that curvature is predicted to always be positive, consistent with previous studies based on mammalian respiration [19], but inconsistent with studies based on plant respiration [16]. Noting the form of $x$, we can identify that $x\propto N$, thus curvature is proportional to $1/N$ and decreases with an increasing number of generations. This network size-based suppression of curvature has been reported in relation to plant respiration [16].

These results demonstrate a current knowledge gap within the field of metabolic scaling, and motivate the simultaneously conducting respiration-based measurements of metabolic rates and vascular imaging to better synchronize prediction with observation.

Appendix F Kernel Density Estimation Results

In this section we present the full results from the kernel density estimation procedure spanning all combinations classification features related to geometric scaling variables.