Secondary osteons scale allometrically in mammalian humerus and femur

Intra-cortical bone remodelling is a cell-driven process that replaces existing bone tissue with new bone tissue in the bone cortex, leaving behind histological features called secondary osteons. While the scaling of bone dimensions on a macroscopic scale is well known, less is known about how the spatial dimensions of secondary osteons vary in relation to the adult body size of the species. We measured the cross-sectional area of individual intact secondary osteons and their central Haversian canals in transverse sections from 40 stylopodal bones of 39 mammalian species (body mass 0.3–21 000 kg). Scaling analysis of our data shows that mean osteonal resorption area (negative allometry, exponent 0.23,R2 0.54,p<0.005) and Haversian canal area (negative allometry, exponent 0.31,R2 0.45,p<0.005) are significantly related to body mass, independent of phylogeny. This study is the most comprehensive of its kind to date, and allows us to describe overall trends in the scaling behaviour of secondary osteon dimensions, supporting the inference that the osteonal resorption area may be limited by the need to avoid fracture in smaller mammalian species, but the need to maintain osteocyte viability in larger mammalian species.


Introduction
The cell-driven process of replacing existing cortical bone tissue with new packets of osteoid (a collagenous matrix which later mineralizes, becoming bone) is termed intracortical bone remodelling and involves the coordinated action of osteoclasts (bone-resorbing cells) and osteoblasts ( In this study, we address this shortcoming. We seek to establish the parameters of any scaling relationships between the dimensions of the intra-cortical bone remodelling process and animal size. To achieve this, we first generated composite images of transverse histological sections of mammalian femora and humeri from a museum collection. We then manually digitized the borders of all intact secondary osteons and their Haversian canals in these images using imaging software. Finally, we measured their areas and analysed these measurements using the standardized reduced major axis method [31] and phylogenetically independent contrasts [32].

Specimens
We obtained access to the historical Bd series of John Thomas Quekett's collection of microscope slides, which today is part of museum collections at the Royal College of Surgeons (RCS) [41]. Quekett's collection is one of the earliest museum collections for comparative histology of bone. Sections of larger species covered only part of the whole bone cross section (figure 1). We selected transverse sections from femur or humerus of mammalian species (total 61 specimens, 56 species; see table 1 (humerus) and table 2 (femur) for a list of specimens with searchable identifiers). The rightmost columns specify whether any secondary osteons were found. Our further analyses dealt solely with specimens that had secondary osteons. For each species, we also list an estimate for adult body mass, obtained from the PanTHERIA database [42] for extant species and from the study by Larramendi [43] for fossil proboscideans, which we used as an input for our scaling analyses. We chose stylopodal bones, as we expected the least variation in loading patterns among species (although this variation may still be considerable) and the greatest number of secondary osteons in these bones [20].
In some cases, the species was unclear, due to missing or antiquated binomial or common names in the historical documentation of the slides [41]. When this occurred, we usually deduced the taxonomy from the description in the catalogue entry. For some species, we additionally resorted to secondary literature (See footnotes of tables 1 and 2) to establish a species with more confidence. Out of the six more markedly uncertain cases, two (Bd 137, 'fossil elephant' and Bd 139, Mastodon sp.) had secondary osteons and were therefore relevant to our main analysis. Analysis of these specimens, as they were fossils, was additionally complicated by possible diagenetic effects. However, scaling parameter estimates calculated excluding these specimens remained comfortably within the overall confidence intervals throughout. Thus, excluding these specimens would not have altered our conclusions and we included them for their comparative value.

Historical details on specimen origin and preparation
The specimens we imaged are likely to have been prepared in different ways by several people. Some display commercial coverslip paper of the time and may have been purchased by Quekett, while others are likely to have been produced in-house by Quekett or his collaborators. In his treatise on the use of the microscope, Quekett describes dry and wet preparation methods for thin sections of bone, and examples of both are present in the collection today [58]. Exact anatomical location of our samples is unknown, although most seem to be taken along the diaphysis, based on seeing smooth endosteal surfaces with no trabeculae, except, as expected, in the two sloth species [59] and some of the marine mammals [60]. Specimens from larger species consisted of only partial cross sections (e.g. figure 1a). Therefore, circumferential position was also unknown.
The origin of the samples used in this study is unclear. Some are likely to have originated from the Zoological Society of London (ZSL), London Zoo, as Richard Owen, for whom Quekett was initially an assistant before becoming his successor, had the right to claim any freshly deceased animal there. Others may have come from skeletons donated to the RCS by John Gould and other nineteenthcentury travelling naturalists. Such donations are described by Quekett in his diaries (e.g. entries of 1 October 1840 and 1 December 1842, accessed at RCS library (MS0027/1)), and sections may have been prepared for the microscope later. Another possible source of specimens is the ZSL Osteological Museum, suggested by a letter from 1840 by ZSL honorary secretary William Ogilby offering this collection to the Hunterian Museum (RCS-MUS/11/1/14, item 52).

Imaging
We imaged the specimens in reflected light fluorescence microscopy using a Leitz Laborlux 12 microscope with a built-in LED light (Excitation 450-490 nm). We manually took overlapping 16-bit grey-scale images, each 1024 by 1024 pixels, covering the entire specimen (4× magnification, 0.12 numerical aperture, pixel size 1.64 µm) with an Orca FLASH4.0LT C11440 digital camera (Hamamatsu Photonics, Hamamatsu City, Japan) and F56-L019 filterset (AHF Analysentechnik, Tübingen, Germany). Consequently, the contrast in the images was generated by green autofluorescence induced by a blue light. In total, we took 3416 images, resulting in an average of 58 image tiles/specimen (minimum 6 (specimen number Bd 91), maximum 180 (Bd 123)). We composed images of entire specimens by applying the Grid/Collections stitching plug-in ( [61], v. 1.2) of the open-source image processing software Fiji ( [62], v. 1.51) to the image tiles.
We found intact secondary osteons in 40 specimens (14 humeri, 26 femora) from 39 species. Using a Wacom digitizer (Wacom, Saitama, Japan), we manually traced the boundary of intact secondary osteons and their Haversian canals in the images. To identify secondary osteons, we followed a method used in previous studies of osteonal dimensions [14,33]. Specifically, this involved classifying a feature as an intact secondary osteon if more than 90% of the cement line and the entire canal boundary were visible, and the aspect ratio (maximum diameter/minimum diameter) was less than 2, to prevent atypical osteon variants from confounding our area measurements. For the same reason, drifting osteons [63] were included if we could draw a trace around the osteon from the non-drifting side of the cement line along a clear lamellar boundary, and excluded otherwise, as previously done in [14]. We found a fracture callus in specimen Bd 270, which we excluded from analysis. We stored the osteon area (the area within the cement line, which is the cross section of the cement sheath), the canal area (the area enclosing the Haversian canal) and the infill area (the difference between osteon area and canal area) for each osteon using a custom ImageJ macro (figure 2). From these directly measured quantities, we also computed two derived quantities for each intact secondary osteon: the infill distance (the distance between the canal border and the cement line, assuming perfectly circular osteons of equivalent area with concentric  Table 1. For each humeral specimen: RCS reference number (searchable on the RCS online catalogue SurgiCat), Quekett's [41] and current terminology, adult body mass estimate [42,43] and whether there were any secondary osteons in it. Secondary literature supporting our assumptions for the terminology can be found in the footnotes.   Haversian canals of equivalent area), and the infill ratio (infill area divided by its corresponding osteon area) using custom ImageJ macros. The number of secondary osteons present in each specimen varied from 1 to 389 (last columns in tables 1 and 2). To prevent our results from being biased by specimens with low numbers of osteons (often small species, some artiodactyls or marine mammals), we repeated our analysis excluding specimens with fewer than 10 measured osteons. This offered some additional robustness to our analyses, because any false positives (features wrongly classified as secondary) would have a particularly high weight in specimens with few osteons. Again, the results for scaling exponents and elevations were within the bounds obtained including all specimens.

Statistical analyses
Scaling analyses aim to describe the relationship between a quantity of interest and body size, usually mass. The quantities of interest for this study included the mean, minimum and maximum of the direct (osteon area, canal area, infill area) and derived (infill ratio, infill distance) measurements described in the previous subsection. As our quantities of interest could not always be assumed to be normally distributed, we distinguished between cases where this assumption was justified and cases where it was not. This was assessed using a Shapiro-Wilk test (with p < 0.05) in R ( [64], v. 3.3.2). We only explicitly state cases where the assumption of normality was unjustified in the results.
To describe the statistical significance of the association between any of our variables of interest and body mass, we report p-values. In all our analyses, we used p < 0.005 and p < 0.05 as cut-off values for high and mild statistical significance, respectively.
To describe the strength of association between any of our variables of interest and body mass, we reported squared correlation coefficients. The squared correlation coefficient was Pearson's R 2 if the data could be assumed to be normal, and the Spearman's ρ 2 otherwise. We interpreted significant relationships with squared correlation coefficients greater than 0.3 as strong, and otherwise as weak.
In cases where the data were not normally distributed (infill ratio), we report the median instead of the mean as a measure of central tendency. In cases where the data were normally distributed, but the fact that osteons which had not finished infilling may have biased our measurements (canal area, infill area and infill distance), we additionally report the median as a measure of central tendency that is robust to outliers.
We assumed our quantities of interest Y obeyed a power law Y = aM b , where M denotes body mass, b the scaling exponent and a the elevation. We obtained estimates, lower and upper bounds for both the scaling exponent b (i.e. the slope on a logarithmic plot) and the elevation a using the standardized (reduced) major axis from the R package SMATR 3 [65] on the base-10 logarithms of our variables of interest. The standardized major axis is the preferred regression method for scaling studies, as it takes into account measurement errors of differing variance on the x-and y-axes [31], and allows testing of hypothesized slopes. We used the same software to test for the null hypotheses of isometry and no correlation. Isometry entails proportional increase with body mass [66], and therefore slope b = 1 3 for the infill distance, b = 0 for the infill ratio and b = 2 3 for all other (area) measurements. Additionally, in order to control for phylogenetic relationships of varying degree between species, we calculated phylogenetically independent contrasts [32] using the pic function of the R package ape species-level phylogenetic tree, including branch lengths [68]. Entries for the extinct proboscideans in our sample were added with a custom R script, using estimated branch lengths from a previous study [69]. SMATR was used to estimate the scaling exponent of independent contrasts, while enforcing zero elevation, as justified in [31]. The rationale behind using phylogenetically independent contrasts is that if the scaling relationships of the contrasts differ strongly from the ordinary scaling relationship for a quantity of interest, phylogeny is an important factor in the study and the data are not as independent as otherwise assumed.
If scaling exponent and elevation estimates calculated for femora and humeri individually were not comfortably within 95% CI of the estimates of combined stylopod data, we report the data separately. Otherwise, we pooled the data for analysis.
Our R code [70] and ImageJ macros [71] are available online.

Direct measurements
Mean osteon area, canal area and infill area showed negative allometric (i.e. less than isometric) scaling (b = 0.23, b = 0.31 and b = 0.22, respectively), as visualized in figure 3. These relationships were different to both null hypotheses (scaling-invariance and isometry) with high significance (p < 0.005). Sample squared Pearson's correlation coefficients (R 2 ) showed that our fits for the means of directly measured variables, where highly significant, explained between 45% and 52% of the variation in our data. The maxima and medians for these quantities showed the same gross scaling pattern (scaling exponents within 0.02 of each other) as the means, with similar values for p and R 2 . The maximum canal area scaled close to isometry, but was still negatively allometric with mild significance (p < 0.05). The variation in our data for minimum areas was considerably less well explained by our fits. The correlation of minimum osteon area and minimum infill area with body mass was weak (R 2 = 0.14 and R 2 = 0.12, respectively, p < 0.05), and was lost completely when excluding samples that had fewer than 10 osteons (p = 0.37). The minimum canal area was not normally distributed (p < 0.005, Shapiro-Wilk), and a Spearman's rank correlation test with body mass yielded a similarly weak relationship (ρ 2 = 0.21, p < 0.005; figure 4). The elevation for the mean canal area was 390 µm 2 , two orders of magnitude less than its analogues for osteon area and infill area, which corresponds to a mean canal diameter of 11.1 µm (assuming circular canals). Estimates and 95% CI for a and b of directly measured quantities are summarized in table 3.

Derived measurements
Shapiro-Wilk tests indicated that the logarithms of the mean, maximum and median infill ratio were not normally distributed (p < 0.005). We therefore performed non-parametric (Spearman's rank) correlation tests on these data, which showed that the mean, median and maximum infill ratio were not correlated with animal body mass (p > 0.18 in all cases). The minimum infill ratio was not normally distributed (p < 0.005, Shapiro-Wilk test), and a non-parametric test showed a very weak correlation with body mass (p < 0.05, ρ 2 = 0.12). The significance of this correlation was lost when analysing the data within separate groups of femoral and humeral specimens. The median infill ratio ranged from 84% to 99% (figure 5). The mean infill distance scaled with negative allometry (b = 0.11, R 2 = 0.44, p < 0.005), indicating that while the infill distance increases on average with species size in absolute terms, larger mammalian species have a relatively smaller mean infill distance (table 4). We observed similar relationships for the median and the maximum infill distance. However, the minimum infill distance is independent of species size, possibly due to the fact we included osteons at an early infilling stage. The overall maximum infill distance (found in specimen Bd 137, straight-tusked elephant) was 180 µm.
These results are summarized in table 4.

Phylogenetic signal
Scaling exponent estimates for phylogenetically independent contrasts of all variables of interest were slightly higher, but in all cases within the 95% CI obtained without any correction for phylogeny (

Discussion
Our analysis shows that cross-sectional parameters of secondary osteons and their Haversian canals scale with negative allometry in the mammalian humerus and femur. Therefore, our key finding is that, where intra-cortical remodelling is present, larger mammals have larger osteonal and Haversian canal areas compared to smaller mammalian species in absolute terms, but smaller osteonal and Haversian canal  Table 3. Summary of the scaling analysis for the osteon area, canal area and infill area: slope estimate with lower and upper bounds (b, b − , b + ), elevation estimate with lower and upper bounds (a, a − , a + ), sample squared correlation coefficient (R 2 ), and p-values for no correlation (p uncorrelated ) and isometric scaling (p isometry ) null.       areas in relative terms ( figure 6 and table 3). The distance between the canal border and the cement line also showed a negatively allometric relationship with animal size, whereas the ratio of the infill area to the osteon area does not scale (table 4). The range of species, clades and adult body mass considered here comfortably exceeds the only previous scaling study concerning itself with secondary osteons [28].
In general, our data support the accepted view that extensive remodelling [20] is present only in larger species. We did not find any secondary osteons in some medium-sized terrestrial species (

Physiological and mechanical constraints on osteon dimensions
Previous studies showed a tight coupling between bone resorption and formation in intra-cortical bone remodelling. Qiu et al. [72] demonstrated a strong (r = 0.97) quadratic relationship between osteon perimeter and infill area in young and elderly humans. Another study [73] showed strong correlations (r > 0.95) between osteon area and infill area for osteons of varying sizes from several species. The size invariance of the infill ratio, and similar exponent estimates for the osteon area and the infill area (differing by less than 0.01) reported here extend this finding to an unprecedented range of mammalian species with secondary osteons. Thus, seemingly, a similar percentage of resorbed bone is filled back in, independent of body size, suggesting that it is crucial for all mammalian limb bones to keep their overall porosity at a similar level, because porosity is an important determinant of bone strength [74].
Sample squared correlation coefficients are lower for the minimum than for the maximum and mean osteon, canal and infill area. Erythrocyte dimensions may form a lower limit for the cross-sectional area of remote blood vessels such as Haversian canals and do not vary with species size in mammals [75]. Idealizing the minimum canal areas as circles resulted in diameters of approximately 10 µm, which are close to the diameter of an erythrocyte (approx. 8 µm, [76]), supporting this notion. Therefore, occasionally, osteons of large species are narrow relative to the entire, inter-species range we report in the present study (and erythrocytes still fit through the capillaries), but wide osteons do not occur in small species.
On the other hand, upper limits for osteon and canal area may differ for species of varying size: The maximum distance from a blood vessel at which osteocytes remain viable has been reported to be 230 µm [16]. This is close to the maximum infill distance we observed in the larger species of our samples and suggests that maintaining osteocyte viability, combined with the previously mentioned size-invariant infill ratio, constitutes a limiting factor for osteon area in large species. For smaller species, however, this limit is never reached, suggesting that a different mechanism dominates the size of osteons in these animals.
The cortical thickness of the specimen with the smallest estimated body mass and evidence of osteons in our sample (red slender loris, Loris tardigradus) ranges between 600 and 850 µm, which is of comparable scale to the estimated diameter of the largest osteon we found overall (580 µm). A resorption area of that size in as thin a cortex is likely to have deleterious consequences for the bone [77]. Creating as large osteons as possible, if an organism can afford the temporary increase in porosity due to a large resorption cavity, may be advantageous from a mechanical point of view after infilling has completed. Larger (and more numerous) osteons have been associated with increased toughness [78][79][80][81]. It has more recently been suggested that larger osteons [82] are geometrically more advantageous for resisting bending and Table 4. Summary of the scaling analysis for the infill ratio and the infill distance: slope estimate with lower and upper bounds (b, b − , b + ), elevation estimate with lower and upper bounds (a, a − , a + ); for the infill distance, sample squared correlation coefficient, and p-values for no correlation (p uncorrelated ) and isometric scaling (p isometry ) null hypotheses are reported. We report the median instead of the mean, and Spearman's ρ 2 , for the infill ratio, as the data were not normally distributed, as indicated by the p-values (p shapiro-wilk ) of normality tests. Isometry was not tested for the infill ratio, as it did not scale in the first place.  compression loads. Studies treating bone as a quasi-brittle material suggest that bone strength and toughness improve with increases in the ratio between dominant inhomogeneity size (i.e. in our case, the typical size of microscopic features such as osteons) and the structural scale (i.e. the bone organ size) [83][84][85]. Therefore, larger osteons, to the extent to which they do not preclude osteocyte viability, in large species may be explained as tissue-level geometrical adaptations to offset mechanical disadvantages that come with increased size, conceptually similar to adaptations at the organ level described in previous studies such as decreasing bending moments with more upright postures [86], increased limb bone robustness [87], cortical thickness and infilling with trabeculae [88], and adapting number and thickness of trabeculae [21].
It has previously been hypothesized that the energy efficiency of the fluid transport to the outermost osteocytes limits the size of secondary osteons [89]. However, that study based its argument on an idealized osteon model, where the velocity of the fluid flow is constant from lamellar boundary to lamellar boundary. A more recent study, using three-dimensional stacks of electron microscopic images and advanced modelling techniques, has shown a highly heterogeneous flow profile in the canaliculi, casting considerable doubt on this assumption [90]. Furthermore, the upper limit to osteon size of 'five to six lamellae', apparently corresponding to 150 µm, proposed in [89] is directly contradicted by observations of up to 13 lamellae in human secondary osteons (e.g. fig. 6 in [91]).
The fact that phylogeny had little effect on our findings is unsurprising, given the diversity of species, clades and ecologies that we sampled from. This further reinforces the notion that our results reflect overall bio-physical limits within which the intra-cortical bone remodelling process operates, and not clade-or species-specific adaptations.

Limitations
The high variability of our data reflects the multitude of species-and specimen-specific characteristics that affect bone remodelling, including age, sex, habitual loading mode, locomotor style, life history, levels of activity, exact specimen body mass, onset of secondary remodelling during ontogeny, environmental influences and more. Specimen age was unknown, but we expect most specimens to be of young adult age, as life expectancy in nineteenth-century zoos was low ( [92], as cited in [93]). There is some debate about the relationship between age and osteon size in humans: while most studies support an inverse relationship (narrower osteons in older people) [94,95], some studies show no effect. Osteon size was not correlated with age in macaques [35]. The effect size of age in inter-species comparisons of secondary osteon dimensions (if there is one) is, however, likely to be smaller in magnitude than the relationship of secondary osteons with mass that we show here. Similarly, we could not take sex into account in our analysis, because it was unknown for all our specimens. Habitual loading directions (compression or tension) may affect osteon size, but it is unclear how. Osteons in the artiodactyl and equine calcaneus and the chimpanzee femur have been shown to be narrower in the compression cortex compared to the tension cortex, but the opposite was true of the equine radius and the human femur [12].
Given that it was difficult to establish the anatomical location and orientation from which the samples were taken, we could not control for this fact in this study. This line of thought is, however, promising for future scaling studies including specimens for which the sample location is known, although caution should be exercised when estimating the loading environment from geometrical parameters [96]. The fact that we sampled only complete osteons may have biased some of the data presented here towards secondary osteons that had not finished infilling. This bias may have been alleviated by sampling Haversian canals from secondary osteon fragments, which are more likely to have finished infilling. However, as some of our samples contained only one generation of secondary osteons (and therefore no fragments), this approach would have severely compromised our ability to compare data throughout species. Our method which captures osteons that had not completed infilling (and ones that stopped infilling at an earlier stage than others [97]), gives a more complete picture relevant to the biomechanical situation. Furthermore, median measurements (which are robust to outliers, including those resulting from incompletely infilled osteons) of the canal area, infill area and infill distance displayed essentially the same scaling patterns as the means throughout, so we are confident that this bias did not affect our conclusions.
Because stylopodal dimensions also scale with body size [24], it seems reasonable to speculate that bone size, and not body size, determines osteon dimensions [98]. However, data for calcanei and femora from equine samples, and metapodials and femora from American black bears, suggest that secondary osteon dimensions are similar across different anatomical sites and across bones of varying size within the same species [73]. This makes it appear more likely that osteon dimensions relate to adult body mass, as our data show, and not to bone size. To distinguish the contributions to osteon size of bone organ size and whole body mass in more detail, a further study, similar to ours, comparing osteon dimensions in bones of the same size across species, and in bones of different sizes within species, would be required.
More insights would be gained if we were to measure the secondary remodelling process in three dimensions. However, while we are now able to image resorption cavities effectively using X-ray microtomography [99], for now, the characterization of cement sheaths in three dimensions involves considerable manual effort [4] and could not have been completed within the time scale of this study.
Using a historical collection of slides presented some further limitations. Owing to the large variation in sample cross-sectional area, unknown anatomical location and unclear boundaries of remodelled areas in some specimens, we did not feel comfortable reporting any osteon areal density (osteons mm −2 ) measurements, as they would not have been comparable between species. Having only a few (one or two) specimens per species further limits our study, as these specimens may not be representative of the species. This could, for example, have been the case if the animal in question suffered from illnesses affecting bone metabolism or regular locomotor activity. To an extent, knowledge of exact anatomical location, precise taxon as well as specimen age, health and provenance was sacrificed in order to sample a large species and size range. We feel this was important, as we would have been prone to concluding that there was no scaling relationship if we had used a smaller number of different species that are easily accessible, given that (for example) humans have relatively wide, and cattle have relatively narrow, osteons for their size. The use of a historical collection additionally avoided the time, effort and difficulties required to obtain and generate an equivalent collection to a modern standard.
Fluorochrome labelling could have helped validate our results. However, this would have constituted only a partial validation: it would have confirmed whether we characterized the secondary osteons formed during the labelling period correctly (i.e. true positives), but would probably have produced a number of false negatives, namely intact secondary osteons formed prior to fluorochrome injection. Using solely labelled osteons would have reduced the number of osteons sampled, increased the bias towards osteons which had not completed infilling at the time of death or had possibly differing infilling rates, and could have therefore confounded our results. Moreover, obtaining labelled material from a large enough range of species would have posed considerable difficulties.

Conclusion
In summary, the images of historically important specimens we produced are available (see data accessibility statement below) to researchers for use in future comparative histological studies. The measurements of secondary osteon and Haversian canal area we presented here have the highest number of mammalian species with the most diversity in clades and the broadest range in size to date. Having such a comprehensive sample was crucial for our analysis, but also intrinsically limited by the uncertainty surrounding historical specimens. We performed a scaling analysis of quantities describing the intra-cortical remodelling process across species, which showed negative allometry for osteon, canal and infill area as well as infill distance, while the infill ratio was size-independent. Our results indicate that minimum osteon dimensions correlate only very weakly with animal body mass and might relate to erythrocyte dimensions. By contrast, the upper limits to osteon dimensions display a strong relationship with adult body size: osteons in small species may be restricted to sizes avoiding fatal temporary stress increases around too large resorption cavities, while osteon size in large species may be dictated by the ability to maintain osteocyte viability.
Ethics. As the samples used were over 150-year-old museum specimens, no ethical approval was required for the present study.
Data accessibility. Full-resolution images can be obtained from RCS for research purposes upon request (low-resolution versions can be found on the RCS online catalogue SurgiCat by searching for the reference numbers listed in tables 1 and 2). Segmentations of the images, and the raw data of the histomorphometrical measurements are available from the Dryad Digital Repository https://doi.org/10.5061/dryad.637r1 [100]. The code used to obtain our results is available on GitHub [70,71].