Uncovering of intraspecies macular heterogeneity in cynomolgus monkeys using hybrid machine learning optical coherence tomography image segmentation

The fovea is a depression in the center of the macula and is the site of the highest visual acuity. Optical coherence tomography (OCT) has contributed considerably in elucidating the pathologic changes in the fovea and is now being considered as an accompanying imaging method in drug development, such as antivascular endothelial growth factor and its safety profiling. Because animal numbers are limited in preclinical studies and automatized image evaluation tools have not yet been routinely employed, essential reference data describing the morphologic variations in macular thickness in laboratory cynomolgus monkeys are sparse to nonexistent. A hybrid machine learning algorithm was applied for automated OCT image processing and measurements of central retina thickness and surface area values. Morphological variations and the effects of sex and geographical origin were determined. Based on our findings, the fovea parameters are specific to the geographic origin. Despite morphological similarities among cynomolgus monkeys, considerable variations in the foveolar contour, even within the same species but from different geographic origins, were found. The results of the reference database show that not only the entire retinal thickness, but also the macular subfields, should be considered when designing preclinical studies and in the interpretation of foveal data.

. Boxplots showing the retina thickness values of the right (a) and left (b) eyes. For each thickness coefficient, numerical data of Mauritius male, Mauritius female, Asian male, and Asian female are plotted. Rectangular boxes represent interquartile ranges (IQR), which extend from Q1 to Q3. A black line in the middle of an IQR indicates a median. Upper whiskers extend to the last datum, which is smaller than Q3 + 1.5 × IQR. Lower whiskers extend to the first datum, which is greater than Q1 − 1.5 × IQR. Data beyond whiskers are outliers and plotted as black circles.
Scientific Reports | (2021) 11:20647 | https://doi.org/10.1038/s41598-021-99704-z www.nature.com/scientificreports/ origin does not seem to play a role, but females tend to have smaller retina thicknesses than males. The data at umbo strongly resemble data at T5. For the retina areas, the patterns are largely the same (Fig. 2). Tables 1, 2, and 3 present detailed summary statistics of all the coefficients. Summary statistics and boxplots are shown for the left and right eyes separately to avoid correlations between the eyes. Table 4 reveal that the retina thickness and corresponding area coefficients (e.g., T1 and A1) are highly correlated (correlations between 0.89 and 0.96). Moreover, columns 9 to 12 in Table 4 show that all four umbo coefficients are highly correlated with T5 (correlations between 0.92 and 0.96). Regarding the correlations among the nine thickness coefficients of T1-T9, Table 5 demonstrates the generally rather high correlation between adjacent coefficients (0.55-0.91). However, for nonadjacent coefficients, this correlation is smaller (0.07-0.85). For statistical analyses, we considered only the nine thickness coefficients T1-T9, excluding the eighth area (A1-A8) and four umbo (TUn, TUt, AUn, and AUt) coefficients.

Statistical analyses. Correlation analysis. Columns 1 to 8 in
Principal component analysis. A principal component analysis (PCA) was used to investigate the latent factors that could explain the variability in the thickness parameters T1-T9. The plots in Fig. 3 reveal that the first three principal components (PCs) explain 87.5% and 87.3% of the variability for the left and right eyes, respectively. That is, the first three PCs mainly determine the values of T1-T9. Table 6 shows the PCA coefficients of the first three PCs for the left and right eyes. The patterns for the left and right eyes are largely the same. The first PC is a weighted average of the nine thickness parameters, with the weights roughly corresponding to the relative size of each thickness. That is, the first PC is an overall thickness factor. The second PC appears to be considering the thickness values at the edges (negative sign) versus those at the center (positive sign). This suggests that the second PC is a center-vs.-edge factor. The third PC seems to be considering nasal versus temporal thickness parameters (mainly reversed signs for PCA coefficients left and right of T5). That is, the third PC is a nasal-vs.temporal factor. The multivariate analysis of variance (MANOVA) ( Table 7) demonstrates that the independent variables of sex and origin have a significant impact on the dependent variables T1-T9 in both eyes. The effect of origin was stronger than the effect of sex. The interaction terms between sex and origin were not significant and, thus, were removed from the final models. Table 8 shows the results of the nine two-way analysis of variance (ANOVA) analyses. Differences in sex are significant for the outer thickness parameters (T1-T3 and T7-T9), whereas differences in origin are significant for inner thickness values (T5 and T6). Interaction terms between sex and origin are not significant and, thus, have been removed from the final models.

Discussion
The main aim of the current study was to generate reference data on a large number of cynomolgus macaque eyes (M. fascicularis) 26,27 , which, because of their genetic similarity to humans, have been successfully introduced into biomedical research [28][29][30][31][32][33] . Alongside histopathologic examinations 31,34 , OCT has been demonstrated as a useful imaging tool to assess ocular toxicity, such as for ocular inflammation, in preclinical studies 24,33,35,36 . Despite all the successes, OCT remains a fairly new technique in animal research, which may account for the lack of reference data. In addition, the analysis of OCT scans from cynomolgus macaques is time-consuming and associated with undesired deviations because of the often manually-thus relatively arbitrary-performed readings of the values 37,38 . Consequently, interpreting such OCT results in the context of the natural variability of macular thickness is impossible.
To overcome these limitations and offer more in-depth knowledge compared with our previous report 39 , the present study successfully implemented an automatic machine learning application that was supplemented with a classic approach in a hitherto unprecedented number of cynomolgus monkeys. Their retinal thickness at the fovea-or the site of highest visual acuity-was assessed to provide reference values. Despite morphological similarities among the eyes of most primates, the obtained data show that in addition to the known differences in the cone numbers [40][41][42] , corresponding variations can also be found at the structural OCT level. These variations occur even within the same species of different geographic origins, which are commonly used in preclinical toxicology studies 43 . This circumstance was also evident in the subanalysis of the umbo area, that is, the site with the highest density of photoreceptors (cones). The integrity of this zone is essential for the best visual acuity, and its disintegration can be the first sign of various diseases, such as age-related macular degeneration, but also, for example, the formation of a macular hole 21 . Our data suggest that despite being the same species, certain morphological characteristics have evolved differently. It is important to be aware of this fact, particularly if, during drug development, researchers decide to alternate trials between monkeys of different origins. Consequently, Table 1. Summary statistics of retina thickness values with respect to eye side, sex, and origin. Values are presented in µm. OD oculus dexter, OS oculus sinister, Std statistical analysis, U umbo, n nasal, t temporal.

Stats
Sex Origin  T1  T2  T3  T4  TUn  T5  TUt  T6  T7  T8  T9 OD www.nature.com/scientificreports/ both origins cannot necessarily be used interchangeably in each use case. This particular difference in macular morphology would be the most relevant if we relied on quantitative macula measurements. The metabolic differences between Asian and Mauritian macaques have been investigated and described previously 42 . As an example, it has been shown that in Mauritius macaques, the retroviral restriction factor TRIM-Cyp is not present, which is in contrast to a higher expression of TRIMCyp in Indonesian macaques, causing methodological challenges for AIDS research 44 . However, so far, no data have been available for ocular structure.
In agreement with a previous study, there was no significant difference in the central depression (T4, T5, and T6) between the sexes 39 . Interestingly, individual ANOVA analyses of the data show significant differences in nasal and temporal thickness values (T1-T3 and T7-T9, respectively) between male and female individuals. This suggests that the averaged values of macula measurements, which are commonly used in ophthalmic research and obtained over a larger range, underestimate local variations, so caution is advised in interpreting these values. Overall, females showed a slightly thinner retina than males-which is similar to humans-although major differences between the provenances and age might also have to be taken into account [45][46][47] .
The mean overall retina thickness in the present study was 199 µm, compared with 244 µm in cynomolgus macaques 39 and 305 µm in humans 48 . This deviation can be explained by the fact that in comparable cynomolgus monkey studies, the value was determined over a larger subfield, which is in contrast to the current subfield study that only calculated the values from the deepest point within the foveola. Regarding minimal foveal thickness, an unpublished study (Vilupuru A. Optical coherence tomography in normal eyes of rhesus monkeys. American Academy of Optometry, 2005) using Stratus time-domain OCT (TD-OCT) measured a value that was approximately 60 µm lower. However, considerable deviations between generations of OCT devices and better reproducibility for the spectral-domain system used in the present study have been reported 49,50 . Furthermore, the number of animals was several times higher in the present study. In summary, a final appreciation of the mentioned TD-OCT data is not feasible because the methodology and values have not been made publicly available in full detail. The current study may be limited because the axial length and refractive status were not known, and diurnal variations have not been measured yet. However, such limits are often inherent in retrospective studies and can be complemented in future prospective studies. Of course, the definition of the reference position nulla can be a point of discussion and may represent an additional limit to naturally occurring variations 51 . Using the suggested fully automatic image processing method, human selection bias with regard to the position of nulla can be largely eliminated or even prevented. Another limitation could be that in the accumulated database, no repeated   Var1  T1  T2  T3  T4  T6  T7  T8  T9  T5  T5  T5  T5   Var2  A1  A2  A3  A4  A5  A6  A7  A8  TUn  TUt Table 5. Pearson correlation matrix for thickness coefficients T1-T9. Adjacent coefficients (e.g., T1 and T2) generally show a higher correlation than nonadjacent coefficients (e.g., T1 and T3). The lower triangle of the matrix has been left empty because the correlation matrix is symmetric. T1  T2  T3  T4  T5  T6  T7  T8  T9   T1  1  www.nature.com/scientificreports/ retina measurements were available. However, for repeated retina OCT measurements, only minor variations were found that were within the resolution and repeatability range of the OCT devices 52,53 . The current number of B-scans may not have captured the effective fovea center in every single eye. In addition, the analyzed region does not correspond to the whole retina, so conclusions can only be derived for the central retina. Nevertheless, the reference database covers the area of greatest cell density 54,55 .
In summary, for the first time, the current study succeeded in applying a hybrid machine learning approach to a large number of macaque eyes to determine the reference values for retinal thickness at the fovea. In addition, the present study has described the morphological variability of retinal thickness and how it relates to sex and the origin of these monkeys. Most importantly, the present study has shown that a thorough awareness of the constraining elements in model species supports the careful selection of the appropriate models for ophthalmic research and appropriate reading of the obtained data. The data provided are important for earlier and more sensitive detection, quantification, and characterization of toxic ocular effects in preclinical safety studies 56,57 . In particular, noninvasive OCT examinations can constitute an additional imaging method for comparative studies between OCT and histopathology 58,59 . However, compared with histology 60 , which is usually performed only once, OCT can be performed within seconds and as often as needed, all without damaging the ocular tissue or bloodstream. Consequently, OCT appears to be an ideal instrument for longitudinal investigations and can presumably enable better characterization and monitoring of lesions. Table 6. Principal component coefficients of the first three principal components for left and right eyes. Principal component analysis (PCA) coefficients of the first three principal components (PCs). The table shows the PCA coefficients for the left (three top rows) and right (three bottom rows) eyes. The patterns are largely the same for the left and right eyes. The first PC is an overall thickness factor with PCA coefficients of T1-T9 corresponding roughly to the relative size of the respective thickness. The second PC is a center-vs.-edge factor. The third PC is a nasal-vs.-temporal factor. The sign of the third PC is mostly reversed for the left and right eyes. Eye  T1  T2  T3  T4  T5  T6  T7  T8  T9 A3218-01)). The protocols for the original drug development studies were reviewed and approved by the respective Institutional Animal Care and Use Committees of their respective contract research organizations. The animals were group-housed in stainless steel cages according to European housing standards (Annex III of Directive 2010/63/EU). The temperature of the animal room was maintained between 20 and 26 °C, with humidity between 30 and 70%. The light cycle was 12 h light and 12 h dark, except during designated procedures. The animals were fed a standard diet of pellets supplemented with fresh fruits and vegetables. Water treated with reverse osmosis and ultraviolet irradiation was freely available to each animal via an automated watering system. Psychological and environmental enrichment was provided to the animals, except during study procedures and activities. The animals were purpose-bred and of Mauritian or mixed Asian origin (the exact geographical location for the latter is unknown).

PC
OCT scans were obtained from 374 eyes of 203 animals. Information regarding age was available for 159 subjects, with an overall mean age of 4.98 years (range 2.5-5 years). Females contributed 147 eyes (39.30%) and males 227 eyes (60.70%); of the total eyes, 186 were left eyes (47.74%), and 188 (50.26%) were right eyes. Regarding the geographic region, 199 eyes originated from Mauritius (53.20%), and 159 eyes were derived from Asia (46.80%). Sixteen eyes were of unknown origin, but the sex was defined.
Inclusion and exclusion criteria. The inclusion criteria were healthy and untreated cynomolgus monkeys derived from Mauritius or Asian origin who were between 30 and 50 months of age and weighed between 2.5 and 5.5 kg. Only healthy eyes, with optically clear ocular media and no observed anterior or posterior segment pathologies, were included in the study. Eyes of an undocumented sex, origin, or eye side were excluded from the subanalyses.

OCT imaging. OCT was performed under general anesthesia with dilated pupils using spectral-domain
OCT on the Spectralis HRA + OCT Heidelberg platform (Heidelberg Engineering, Heidelberg, Germany) as reported 37,39 . The same scan protocol was selected for all animals with a horizontal line scan pattern (centered over the fovea) with a size of 20° × 20°, 25 raster lines separated by 221 μm (scan length 5.3 mm, 512 × 496 pixels, scan depth 1.9 mm). The data were exported directly from the OCT device as original B-scan files in a bitmap image data (BMP) format using the manufacturer's software.
The fovea position-location problem. Because the examined animals for obvious reasons could not follow the instructions of the operator to exactly and steadily fixate on a presented target, the determination of the fovea center by an OCT operator could be relatively arbitrary. The use of subjective criteria for the definition of the fovea may affect the preciseness of a measurement 61 . Although the term fovea is frequently used in the clinical context, it does not represent a precise anatomical term or position, so the determination of fovea landmarks may be associated with uncertainties 55,62 .
Currently, a concaviclivate fovea (fovea from the Latin "ditch" or "pit") represents a depression of the center of the retina 9 (Fig. 4). Thereby, a supposed margin can be assumed, which transits into a descending clivus or slope (clivus from the Latin "slope" or "hill"), resulting in a certain bowl-shaped configuration. The base is referred to as the foveola or "little fovea. " Although the umbo is interpreted clinically as the center of the fovea, such a definition proves to be rather problematic in cross-sectional imaging in animals because the geometry and extension of the fovea differ remarkably among different species 55 . Definition of nulla as an OCT structure-based fovea parameter. The fovea can be visualized very well with an excavation of the central retina on OCT imaging (Fig. 4a). To eliminate the position-location problem of the fovea toward an objective assessment and in a slight customization to previous reports 51, 55 where the deepest point of the foveal pit was indicated, a modified interpretation of the deepest position within the retina is proposed based on structural OCT information.
For this purpose, an automatic fovea depression contour finder approach was pursued, in which the definition of the particular fovea center (displayed on the OCT-B scan) was defined as the deepest point within the foveolar cavity (Fig. 4b). The resulting position was referred to as the nulla, which implies that on the particular B-scan or a stack of B-scans no deeper excavation could be present below the nulla within a fovea. The retinal thickness at the nulla position corresponds to the thinnest retinal thickness within the pit or inner fovea, where the incident light can interact most directly with the photoreceptors. This position holds an important value because the area within 150-200 µm around the nulla represents a crucial site containing the highest concentration of photoreceptors (cones) 55  . The first step is the generation of pixel-wise semantic segmentation maps of the retina compartment using a convolutional neural network (CNN) for each B-scan (Fig. 5). The CNN was developed and described in detail in a previous study 63 . In summary, it uses a modified U-Net architecture 64 with 22 convolutions, 5 transposed convolutions, and 5 skip connections, which have previously been shown to also be highly effective in learning The area of the fovea is characterized as a central depression (single arrows) descending in a more or less symmetrical curvilinear shape (double arrowheads) to the bottom, which is named the foveola. The deepest location inside the fovea is determined as the nulla (marked as a red dot) and represents an OCT-based anatomical landmark. The umbo in the OCT image is designated as the proposed center of the fovea (highlighted as the green area between the two arrows) near the nulla at a distance of 100 µm to each side. Based on these conventions, the following OCT-based sequence can be proposed: macula lutea > perifovea > parafovea > fovea > foveola > umbo > nulla. The retinal pigment epithelium is highlighted in blue.  (2) adding a random rotation between − 8° and 8° degrees, increasing the training set size to 2025 B-scans. Regarding CNN training, the model parameters were initialized 65 and learned by minimizing an unweighted pixel-wise cross-entropy loss summed over the entire input. Adam optimization was used on a single NVIDIA TITAN-X GPU 66,67 . An initial learning rate of 6 × 10 −5 and a mini-batch size of eight images was chosen, which had proven suitable in preliminary experiments. Training was stopped after 1920 iterations (7.6 epochs) after the validation set accuracy reached a plateau. On the test set, the differences between the CNN's predictions and the annotations of the three human graders were, on average, smaller than the intergrader differences. The algorithm's input are B-scans (Fig. 5a), here rescaled to 512 × 512 pixels and with its output corresponding pixel-wise semantic segmentation maps of the compartments vitreous, retina, choroid, and sclera (Fig. 5b). For the current study, only the retina compartment segmentations were further processed (Fig. 5c). The retinal inner boundary was defined as the transition from the hyporeflective vitreous to the hyperreflective retinal nerve fiber layer, more specifically in the area of the ILM. The retinal outer border was defined as the outer demarcation zone of the hyperreflective retinal pigment epithelium, just above the hyporeflective zone of the choriocapillaris. Detailed CNN architecture is shown in Fig. 2  Segmentation artifact removal. The second algorithmic step is the removal of segmentation artifacts.
A small number of semantic segmentation maps contained small patches of misclassified regions, that is, segmentation artifacts. The CNN output in Fig. 6 shows an example of two such misclassified choroid regions in the sclera compartment. These segmentation artifacts were corrected by the following approach: in a segmentation map, the four main compartments (vitreous, retina, choroid, and sclera) are expected to form large, connected regions. This is true even in the presence of segmentation artifacts. Thus, the four largest connected regions generally correspond to the vitreous, retina, choroid, and sclera compartments. If additional connected regions are present, they must represent segmentation artifacts. If a segmentation artifact is completely surrounded by a compartment, the artifact can be removed by replacing it with the label of the surrounding compartment. Because the segmentation artifacts were small and did not occur at the compartment borders, they were effectively removed using this approach.
Nulla identification. The third algorithmic step is the identification of the ILM (the border between vitreous and retina) from the semantic segmentation map (Fig. 6). This approach generally yielded noisy ILMs, which did not allow for a reliable identification of the lowest point: the nulla. Therefore, in a fourth algorithmic step, the identified ILM was smoothed using a moving average, with the moving average window being applied in two dimensions simultaneously: in the B-scan width dimension (window size 11) and B-scan "stack" dimension (window size 5). Afterwards, the fifth algorithmic step identified the lowest point on the smoothed internal limiting membrane (Fig. 6). If multiple lowest points were found (usually adjacent), the coordinates of their center of mass were used as the lowest point and, if necessary, rounded to the coordinates of the voxel nearest to it. By utilizing this approach, the "central" B-scan was identified, along with the pixel coordinates of the nulla. Figure 6. Illustration of the image processing pipeline that consists of six steps. First, a raw B-scan is used as input to a modified U-Net 64 . The output is a semantic segmentation map. In between, 11 network layers are schematically visualized, including the network's five skip connections (indicated by black arrows) 61

Retinal thicknesses and areas.
The sixth algorithmic step is the measurement of retinal thicknesses and its areas. Using the semantic segmentation of the retina compartment on the central B-scan and the position of the anatomically and structured OCT-determined nulla, an imaginary line was orthogonally positioned in relation to the underlying retinal pigment epithelium to determine the axial diameter (Fig. 7). Subsequent measurements of retinal values were made at intervals of 500 µm to the side, up to a maximum of 2000 µm distance from the nulla. This allowed for measurements of nine retinal diameters (marked as thickness parameters T1-T9) in the axial direction, as well as eight intervening retinal areas (A1-A8), hence providing a total of 17 parameters to quantify the retina features (Fig. 7). An additional subanalysis in the important area of the umbo was performed to determine the parameters for the most central receptors, which are responsible for the best visual acuity (central bouquet of cones) 41,51 . For this purpose, the thickness values and retina intervening surface areas were each determined laterally at an interval of 100 µm to the nulla. Thus, four further parameters were added: one additional nasal thickness (TUn), one temporal thickness (TUt), and two additional retinal surface areas from the nasal (AUn) and temporal (AUt) retina areas. Finally, including T5 (also indicated as the nulla), the umbo subfield analysis used a total of five parameters.
The first algorithmic step was performed in Python v3.5 and TensorFlow v1.14.0 69 . Steps two to six were performed in C# (v7.0, NET Framework v4.6). The combination of CNN-based semantic image segmentation with a fully automated fovea-finding algorithm, which is based on classical computer vision, can be described as "hybrid image processing. " Data analysis. Summary statistics and visualization. The summary statistics of the mean, minimum, and maximum were calculated for each of the measured retinal thickness and area coefficients on subsets of the data. Boxplots were used to visualize the distribution of the data and compare the two groups with each other (e.g., male vs. female). For the nulla, its retina diameter (T5), its adjacent retina surfaces (A4 and A5), and the average mean values were calculated for all eyes and both sexes.
Statistical analyses. Statistical analyses were performed (1) to investigate the observed variability in the data and (2) investigate differences with respect to sex and origin. The statistical analyses consisted of four parts.
First, a Pearson correlation analysis was performed to investigate the correlation among the 11 thickness and 10 area coefficients. As a result of the high correlation between (1) the thickness and area coefficients and between (2) the umbo subanalysis coefficients and T5, it was decided to include only the nine thickness coefficients T1-T9 in subsequent statistical analyses. For the Pearson correlation analysis, the left and right eyes were combined. Separate analyses of the left and right eyes yielded results almost identical to the results of the left and right eyes combined.
Second, a principal component analysis (PCA) was performed separately for the left and right eyes 70 . PCA does not provide insights into the effects of the independent variables of sex and origin but can identify the latent factors underlying the variability observed in T1-T9. Scree plots were used to visualize the eigenvalues of each PC. Each eigenvalue corresponds to the amount of variability explained by the corresponding PC. PCA coefficients are shown for the first three PC.
Third, a two-way MANOVA 71 was used to jointly investigate the effects of the independent variables of sex and origin on T1-T9. MANOVA assumptions were checked with diagnostic plots (not shown). One multivariate outlier was removed from the group of right eyes. Wilks' lambda was used to measure the impact of sex and origin. The MANOVA was separately calculated for the left and right eyes. One-and two-dimensional retina parameters are illustrated in a left retina. From the anatomically deepest location within the foveolar cavity (marked as nulla, red dot) and orthogonally to the retinal pigment epithelium, side-by-side measurements were taken at 500 µm intervals starting nasally toward the temporal retina with respect to the axial retina diameters (a) (shown as pink lines and marked as retinal thickness parameters T1-T9) and the retinal surface areas (b) in between (highlighted in red and marked as areas A1-A8). Subfield analyses included the regions at a distance of 100 µm from nulla (highlighted with white lines and umbo surface area in green). The same procedure was performed for all eyes. For a better understanding, the segmented retina is depicted as a blue overlay. Scale bars, 500 µm. www.nature.com/scientificreports/ Fourth, two-way ANOVA tests were used to investigate the effects of sex and origin for each coefficient T1-T9 individually. The coefficients T1-T9 are obviously correlated with each other and, thus, not independent. Consequently, the ANOVA results are not independent of each other, and the p-values might be inaccurate. Nevertheless, we decided to perform individual ANOVA analyses because it allowed us to gain insights into which parts of the retina are responsible for the differences between male and female and Asian and Mauritius monkeys. To adjust for multiple testing, we performed Bonferroni corrections by dividing the significance levels by the number of tests, which was nine in our case. Variables that were significant at p < 0.001/9 are indicated with "***. " Variables that were significant at p < 0.01/9 are indicated with "**, " and variables that were significant at p < 0.05/9 are indicated with "*. " Finally, variables that were significant at p < 0.1/9 are indicated with ". ". ANOVA assumptions were checked with the diagnostic plots (not shown). ANOVA tests were performed separately for the left and right eyes.
The 374 eyes contained 16 eyes of an unknown origin, which were excluded from the MANOVA and ANOVA analyses. Some monkeys contributed a left eye and right eye. Thus, the left and right eyes are not independent of each other. Consequently, PCA, MANOVA, and ANOVA analyses were performed for the left and right eyes separately. All summary statistics, data visualizations, and statistical analyses were performed in Python v3.8.5. PCA and statistical tests were performed with the Python package statsmodels v0.12.1. Visualizations were generated using the Python package Matplotlib v3.3.2.

Data availability
The measurement datasets generated and analyzed during the current study are included in this manuscript and made available as Supplementary Data S1.