Residential Exposure to Urban Trees and Medication Sales for Mood Disorders and Cardiovascular Disease in Brussels, Belgium: An Ecological Study.

BACKGROUND: The available evidence for positive associations between urban trees and human health is mixed, partly because the assessment of exposure to trees is often imprecise because of, for instance, exclusion of trees in private areas and the lack of three-dimensional (3D) exposure indicators (e.g., crown volume).
OBJECTIVES: We aimed to quantify all trees and relevant 3D structural traits in Brussels (Belgium) and to investigate associations between the number of trees, tree traits, and sales of medication commonly prescribed for mood disorders and cardiovascular disease.
METHODS: We developed a workflow to automatically isolate all individual trees from airborne light detection and ranging (LiDAR) data collected in 2012. Trait data were subsequently extracted for 309,757 trees in 604 census tracts. We used the average annual age-standardized rate of medication sales in Brussels for the period 2006 to 2014, calculated from reimbursement information on medication prescribed to adults (19-64 years of age). The medication sales data were provided by sex at the census tract level. Generalized log-linear models were used to investigate associations between the number of trees, the crown volume, tree structural variation, and medication sales. Models were run separately for mood disorder and cardiovascular medication and for men and women. All models were adjusted for indicators of area-level socioeconomic status.
RESULTS: Single-factor models showed that higher stem densities and higher crown volumes are both associated with lower medication sales, but opposing associations emerged in multifactor models. Higher crown volume [an increase by one interquartile range (IQR) of 1.4×104 m³/ha] was associated with 34% lower mood disorder medication sales [women, β=-0.341 (95% CI: -0.379, -0.303); men, β=-0.340 (95% CI: -0.378, -0.303)] and with 21-25% lower cardiovascular medication sales [women, β=-0.214 (95% CI: -0.246, -0.182); men, β=-0.252 (95% CI: -0.285, -0.219)]. Conversely, a higher stem density (an increase by one IQR of 21.8 trees/ha) was associated with 28-32% higher mood disorder medication sales [women, β=0.322 (95% CI: 0.284, 0.361); men, β=0.281 (95% CI: 0.243, 0.319)] and with 20-24% higher cardiovascular medication sales [women, β=0.202 (95% CI: 0.169, 0.236); men, β=0.240 (95% CI: 0.206, 0.273)]. DISCUSSION: We found a trade-off between the number of trees and the crown volumes of those trees for human health benefits in an urban environment. Our results demonstrate that conserving large trees in urban environments may not only support conservation of biodiversity but also human health. https://doi.org/10.1289/EHP9924.

Introduction

Urban green spaces, and especially urban trees, have been associated with numerous human health benefits through a number of regulating and cultural ecosystem services (Aerts et al. 2021; Chen et al. 2019). Exposure to trees has been associated with better mental health (Marselle et al. 2020), improved respiratory health (Rao et al. 2014), lower rates of cardiovascular disease (Donovan et al. 2015; Knobel et al. 2021), lower odds of diabetes and hypertension (Astell-Burt and Feng 2019), lower odds of sleep disorders (Astell-Burt and Feng 2020a), improved health via improved physical activity (Ulmer et al. 2016), and decreased mortality (Barboza et al. 2021).

A recent scoping review of 18 studies on green space interventions and cortisol outcomes demonstrated the potential for a positive association between active and passive exposure to green space and stress reduction (Jones et al. 2021). For instance, in a study in 1,027 children included in the Dutch Prevention and Incidence of Asthma and Mite Allergy (PIAMA) birth cohort, residential exposure to green space in a buffer of 3,000 m, measured as average normalized difference vegetation index (NDVI; an indicator of greenness, computed as the difference between near-infrared and red reflectance divided by their sum) and total relative cover of green space, was associated with lower stress levels as indicated by larger diurnal decreases in saliva cortisol levels (Bloemsma et al. 2021). Psychological, physiological, and environmental stress relief may be the common mechanism behind the beneficial effects of trees and urban green space on various aspects of human health, in particular, stress-mediated conditions, such as cardiovascular disease and mental disorders (Honold et al. 2016; Lanki et al. 2017; Markevych et al. 2017; Nieuwenhuijsen et al. 2017; Zhang et al. 2021).

Trees may reduce heat stress, relieving the urban heat island effect, through leaf evapotranspiration and shading (Speak et al. 2020), and tree crowns may help to reduce stress caused by noise from traffic and other urban activities (Tashakor and Chamani 2021). Some reports suggest that trees with some species-specific traits may help to reduce air pollution by uptake and capture of air pollutants, thus helping reduce oxidative stress (Grote et al. 2016; Roy et al. 2012; Salmond et al. 2016). Therefore, tree crown volume can play a large role in the magnitude of stress relief, given that it indicates the area of ground that can be shaded and the total leaf area that determines the magnitude of leaf evaporation, transpiration, and interception (Kong et al. 2016; Lin and Lin 2010; Pretzsch et al. 2015; Rafiee et al. 2016; Shahidan et al. 2010; Speak et al. 2020; Yin et al. 2011).

Previous studies on the impacts of urban trees on human health have used satellite imagery derived greenness (NDVI), satellite imagery derived land use and land cover (LULC) maps (Aerts et al. 2020; Browning and Rigolon 2018; Dzhambov et al. 2018), inventories of street trees and trees in parks and other public spaces (Kardan et al. 2015; Marselle et al. 2020; Rugel et al. 2017; Taylor et al. 2015), or street photography (Helbich et al. 2019; Li and Ghosh 2018; Liu et al. 2020; Nguyen et al. 2018; Villeneuve et al. 2018) to quantify residential exposure to urban trees. These studies have found mixed effects of urban trees on human health, and this may be related to the varying quality of the exposure assessment. NDVI is not optimal to quantify exposure to trees because it does not differentiate trees from other vegetation efficiently. LULC maps such as the CLC (CORINE Land Cover) (EEA 2019) data set may have been too coarse (1:100,000), with small green spaces and individual trees masked in mixed or nongreen land use pixels (Aerts et al. 2020; Labib et al. 2020). Nevertheless, small green spaces and trees have been demonstrated to have significant associations with human health improvement (Wood et al. 2017), and therefore their omission may affect exposure–health outcome studies. Although more detailed LULC data (e.g., 1:2,000; Hooper et al. 2018) could be used to identify smaller greenspaces, some individual trees might still not be registered. Studies based on inventories of street trees and trees in other public spaces have excluded trees in private spaces (Marselle et al. 2020), and seasonality may have had a substantial impact on greenness derived from street photographs (Pelgrims et al. 2021). As a consequence, exposure to urban trees may have been underestimated or biased (e.g., misclassification between trees and other vegetation types) in these studies. In addition, most of the exposure indicators derived from different data sources were two-dimensional (2D) in their representations (e.g., tree canopy cover), whereas 3D representations (e.g., tree crown volume) remain less explored, possibly leading to inaccurate indication of associations between urban trees and human health (Jiang et al. 2017).

Airborne light detection and ranging (LiDAR) imagery, comprising a dense 3D point cloud with precise height information of the targets, has enabled more accurate isolation of urban trees from grass, shrubs, and other vegetation types and has been used in studies on the impacts of urban trees on human health (Jarvis et al. 2020; Reid et al. 2017; Ulmer et al. 2016). LiDAR also facilitates individual tree crown delineation by using existing segmentation algorithms, something that is critical for automatic inventory of tree density in all public and private areas (Alonzo et al. 2014; Kwong and Fung 2020). Individual tree crown delineation is also the basis for estimation of 3D tree traits such as height, crown base height, and crown volume using remote sensing. These tree traits, together with tree density, can provide a more comprehensive and reliable quantification of exposure to urban trees compared with 2D tree cover data (Zhao et al. 2021), but they have not been explored in previous studies on the impact of urban trees on human health.

Therefore, the first aim of this study was to use high-resolution airborne LiDAR data to improve the assessment of exposure to urban trees and a number of 3D traits related to their potential health effects. Second, we aimed to investigate the associations between reliable residential urban tree exposure and indicators of stress-related disease by analyzing sales data of medication commonly prescribed for cardiovascular disease and mood disorders.

Methods

Study Design and Population

The present study is part of the GRESP-HEALTH (Impact of green/blue spaces on specific morbidity and cause-specific mortality in Belgium) project (Casas et al. 2015), which is a nationwide ecological study on the association between residential exposure to green spaces and specific morbidity in Belgium. The present study focusses on the Brussels Capital Region (161 km2), the administrative region and urban agglomeration encompassing 19 municipalities, including the capital of Belgium, the City of Brussels (hereafter Brussels). Data in this study were analyzed at the census tract level. Census tracts are the official and the smallest administrative spatial units for statistical purposes at higher resolution than the municipality (“statistical sectors”), defined by the Belgian Statistical Office (Statbel). The total number of census tracts in Brussels is 722, with an average census tract surface area of 22.5 ha (; range: 1.4–754.6 ha) and an average of (; range: ) (Statbel 2021).

Medication Sales Data

This study used health care data from the Belgian social security agency InterMutualistisch Agentschap-L’Agence InterMutualiste (IMA-AIM). The IMA-AIM manages health care data collected by the seven Belgian health insurance funds. In Belgium, health insurance is mandatory and the population in the IMA-AIM database corresponds to about 98% of the Belgian population (as registered in the national register). The IMA-AIM provided data on reimbursed medication sales delivered by pharmacies for cardiovascular disease and mood disorders, both yearly, from 2006 to 2014. The data represented the number of adults (19–64 years of age) per census tract and per year for whom at least one refundable medication was prescribed at least once during the study period (“prescription–patients”). The data were provided by sex. General cardiovascular medication was defined as all reimbursed drugs included in the Anatomical Therapeutic Chemical (ATC; https://www.whocc.no/atc) codes B01A (antithrombotic agents), C01 (cardiac therapy, including cardiac glycosides, antiarrhythmics, cardiac stimulants, and vasodilators), C02 (antihypertensives), C03 (diuretics), C07 (beta blocking agents), C08 (calcium channel blockers), and C09 (agents acting on the renin–angiotensin system) (see also Aerts et al. 2020). Medication for mood disorders was defined as all reimbursed drugs included in the ATC codes N05 (psycholectics) and N06 (psychoanaleptics) (see also Aerts et al. 2022). The aggregated data did not contain information on the frequency of use. Census tracts with no more than five reimbursed persons or registered inhabitants in at least 1 y during the study period were excluded by the IMA-AIM because of privacy reasons (Figure 1). All health data were used under the license of the IMA-AIM. The protocol for this study did not require ethics approval or consent to participate.

Figure 1.

Flow diagram of the census tracts in Brussels, Belgium, included and excluded in the study. Note: ASMR, age-standardized medication sales ratio; , number of census tracts; , minimum number of reimbursements between 2006 and 2014 (at the level of the census tract, unstandardized); POP, population (number of inhabitants).

Exposure to Urban Trees

Lidar data collection and preprocessing.

Airborne LiDAR data, with an average point density of , were acquired in winter 2012 and provided by the Brussels Regional Informatics Center (CIRB; Degerickx et al. 2018). After removing noise points (i.e., isolated points) using LAStools (version 171215; Rapidlasso GmbH) software, a digital surface model (DSM) with a spatial resolution of was derived from the highest LiDAR point in each pixel. A canopy height model (CHM), further used for individual tree delineation, was generated using the difference between DSM and a digital terrain model (DTM; provided by the CIRB). The CHM was further smoothed (i.e., median filtering) to fill in small gaps. A normalized LiDAR point cloud, further used for estimation of tree structural traits, was obtained by calculating the height difference between the LiDAR point cloud and the DTM. The processing of the LiDAR data was conducted in OPALS software (version 2.4.0; Otepka et al. 2012).

Tree crown identification and delineation.

A workflow that included object-based classification, individual tree delineation, and screening incorrect trees was adapted from Degerickx et al. (2018) and Chi et al. (2020) to obtain individual tree crown polygons in the study area. In the object-based classification, a tree index image with a spatial resolution of was generated from the normalized LiDAR point cloud by calculating the height difference between the first and last return in each pixel (O’Neil-Dunne et al. 2014). This index is based on the fact that LiDAR pulses can penetrate tree crowns but not buildings, enabling the differentiation of tree objects from building objects. Objects with a tree index and a height (identified from the CHM) were initially labeled as trees. Individual tree delineation was then conducted on large tree objects using a watershed segmentation algorithm (Wang et al. 2004). The object-based classification and individual tree delineation was performed in eCognition (version 9.4; http://www.ecognition.com/).

The resulting tree crown polygons contained misclassified trees, such as building edges, that could have a high tree index value. Degerickx et al. (2018) used an NDVI image with a spatial resolution of to remove misclassified trees. However, because of a slight geometric mismatch between their NDVI image and LiDAR data, quite a number of nontree objects were still retained in the resulting individual tree crown data set. In addition to height information, LiDAR also provides intensity information that reflects the material characteristics of objects and which has been used for urban land cover classification (Zhou 2013). To find the threshold for the intensity that can differentiate trees from building edges, we randomly selected 2,000 tree crown polygons and 2,000 building edge polygons based on the CHM and calculated the average intensity value for the LiDAR points within each polygon. We found that 98% of selected tree crowns were retained and 84% of the selected building edges were removed when using an intensity threshold value of 250 digital numbers (DNs). By applying this threshold to all delineated trees in the 722 census tracts, 616,379 (76% of all the delineated trees) correctly delineated trees were identified and used for further analysis. The screening of incorrect trees was conducted in Python (version 3.8.5; Python).

A field-measured tree database, including street trees and trees in parks with the coordinates of tree stems, was used to assess the tree delineation accuracy. We randomly selected 2,000 trees from the database and the proportion of field-measured tree stems contained exactly in one delineated tree crown polygon was calculated.

Derivation of indicators of exposure to trees.

Tree height, crown base height, and crown volume were derived from the normalized LiDAR point cloud for each tree. Tree height (H) was defined as the height of the highest point in the tree point cloud. We developed an index () to estimate crown base height and extract crown points. The index was based on the fact that when a tree point cloud is sliced into bins with a same height (e.g., ) along the z-axis, bins from tree crowns contain more points than those from tree trunks. The index was thereby calculated as a quotient of the number of LiDAR points in a bin and the number of all points from the tree. We again relied on the 2,000 correctly delineated trees (see the section “Tree crown identification and delineation”) to determine the threshold for that separates crown points and noncrown points. First, LiDAR points with a height were extracted for the 2,000 trees. We assumed that these extracted points were from tree crowns. Next, for each of the 2,000 trees, the average for the extracted crown points () was calculated using . was the number of LiDAR points with a height and was the number of all LiDAR points. The threshold for that separates crown points and noncrown points was finally set as the mean for the 2,000 trees (mean ) and applied to all the 616,379 correctly delineated trees. Briefly, starting with the bin with the height of (excluding the lowest bin containing ground points), the height of the first bin with a was regarded as crown base height (). LiDAR points with a were labeled as crown points and used to calculate crown volume using a 3D convex hull algorithm in package scipy in Python (Figure 2). For each tree point cloud, the algorithm identifies the outmost points, which are then triangulated with Delaunay triangulation and meshed for surface generation. The LiDAR data used in the present study had a high point density and were able to record detailed crown structure in winter. Standard deviations (SDs) of H and , indicating the variation in tree structure, were also calculated because recent studies have suggested that human health may be associated with tree (structural) diversity (Aerts et al. 2018; Marselle et al. 2021).

Figure 2.

Estimation of tree crown volume using LiDAR data and 3D convex hull algorithm: (A) LiDAR points for a delineated tree and (B) the reconstructed tree crown surface. 3D convex hull algorithm identifies the outmost points (in black), which are triangulated with Delaunay triangulation and meshed for surface generation. Note: 3D, three dimensional; LiDAR, airborne light detection and ranging.

Potential Confounders

An earlier study in Belgium found that socioeconomic deprivation is a strong predictor for cardiovascular medication sales (Aerts et al. 2020), which was consistent with other studies that have reported that socioeconomic confounders may have a higher impact on (cardiovascular) health than exposure to urban green space (Kabisch 2019). Therefore, we included sociodemographic characteristics as potential confounders in the analyses. Sociodemographic data were provided by Statbel. Data were derived from the 2001 census. We included three indicators of socioeconomic status (SES) that are related to socioeconomic deprivation (as in Aerts et al. 2020): the percentage of foreign-born inhabitants from lower- and mid-income countries (LMICs; defined as those countries that did not belong to the World Bank classification of “high-income economies” in 2000), the percentage of unemployed inhabitants (employment deprivation), and the percentage of primary educated or lower among the 25- to 64-year-old inhabitants (education skills and training deprivation). Census tracts with missing SES data (i.e., not provided by Statbel) were not included in the analyses (91 census tracts) (Figure 1).

Statistical Analyses

Age-standardization of the medication sales data was essential to produce comparable measures of prescription behavior because the number of prescriptions in a given area depends on the age structure of the area (Aerts et al. 2020, 2022). The number of observed prescription–patients in an area was divided by the expected number of prescription–patients in that area if the prescription behavior were the same as that of the overall Belgian population, taking into account the age- and sex-specific structure of the area. This standard medication sales ratio was then adjusted with the crude rate of prescriptions (Aerts et al. 2020, 2022). The age-standardized medication sales were expressed as the number of reimbursed prescription–patients per 1,000 inhabitants (as a positive integer). We used the mean annual rate of medication sales for the period 2006 to 2014 as the outcome variable (the overall mean for exploratory analyses, and for women and men and for both medication types separately for our detailed analyses).

Preliminary analyses demonstrated that the privacy-restricted census tracts ( reimbursed persons in at least 1 y) had considerably higher tree crown volumes () than the initially included census tracts (included tracts: vs. privacy-restricted tracts: ; independent samples -test , ). To avoid the introduction of sampling bias by excluding census tracts with few medication sales and high exposure to urban trees, we assigned 40 prescriptions/1,000 inhabitants per year to privacy-restricted census tracts with inhabitants (120 census tracts with a total of inhabitants). Census tracts without inhabitants (among others, forests, parks, railway infrastructure, cemeteries, and industrial areas without inhabitants) were excluded from the analysis (27 census tracts) (Figure 1). The value assigned to restricted tracts with inhabitants was close to the threshold value for privacy restriction given that the lowest age-standardized medication sales ratio (ASMR) reported in the IMA data of medication for cardiovascular disease and mood disorders was 42 prescriptions/1,000 inhabitants per year. If the true rate of prescriptions would have been higher, these census tracts would not have been privacy restricted. The assigned value was thus an estimation of the maximum rate that could have been observed in the restricted tracts, and this approach is therefore conservative.

Generalized linear models based on the Poisson distribution with a log-link function were used to investigate the associations between exposure to urban trees and medication sales. First, we explored single-factor models with the overall average medication sales (both medication types for men and women pooled) as dependent data and stem density and crown volume separately as independent variables. Then we obtained estimates for a two-factor model (crown volume and tree density in the same model) and for a SES-adjusted model for the overall average medication sales. In our main models, all tree exposure variables (i.e., stem density, crown volume, and structural variation) were included in the same model along with the SES factors, but separately for men and women and medication types. Models were fit for men and women independently because medication sales for men and women at the level of the census tract were separate response variables (see also Aerts et al. 2022). Because the medication sales data were aggregated at the census tract level and there was no information on subject-level medication use, we were not able to study combined medication sales. The number of included census tracts was 604 (Figure 1). Because SES variables have consistently been found to be significant in explaining associations between green space and human health in our study area (Aerts et al. 2020; Pelgrims et al. 2021; Trabelsi et al. 2019), only fully adjusted models (i.e., all SES variables included in the same model) were evaluated for our main models. Tree density and tree crown volume were entered as increments in interquartile ranges (IQRs) to obtain comparable parameter estimates. In addition to the main models, we also calculated estimates stratified by population density and SES factors, using for each stratification variable the median as the threshold for creating the strata. Statistical analyses were performed with IBM SPSS Statistics Subscription (version 11-2018) software.

Results

Census Tract Characteristics

The characteristics of the census tracts that were included in the study are presented in Table 1. Average percentages of foreign-born inhabitants from LMICs (%LMIC), unemployed inhabitants (%Unemployed), and primary-educated-or-lower inhabitants (%Low education) were 16.3, 19.7, and 16.5, respectively.

Table 1

Characteristics of the census tracts in Brussels, Belgium, included in the study () stratified by quartiles of crown volume ().

Categories	Characteristics	All census tracts		Tree crown volume quartile (×10^4 m^3/ha)
		Mean±SD	Range	Mean±SD (Q1: 0–0.47)	Mean±SD (Q2: 0.48–0.89)	Mean±SD (Q3: 0.90–1.63)	Mean±SD (Q4: 1.64–5.65)
SES	%LMIC	16.3±9.8	2.4–48.9	24.7±10.2	17.5±8.8	12.4±6.2	10.5±6.5
	%Unemployed	19.7±10.3	4.4–56.3	28.2±11.0	20.4±8.5	16.2±7.8	14.0±7.5
	%Low education	16.5±11.3	1.0–48.9	26.6±11.6	17.3±9.5	12.8±8.3	9.2±7.2
Treesa	Trees (n)	513±919	5–14,227	134.6±186.2	279.5±164.3	474.5±557.3	1,162.7±1,547.3
	Trees (n)/ha	24±14	0–91	9.1±4.1	17.6±4.3	26.1±5.6	43.0±10.7
	Crown volume (×10^4 m^3/ha)	1.2±1.0	0–5.6	0.3±0.1	0.7±0.1	1.2±0.2	2.5±0.8
	Height variation (SD)	5.5±1.3	1.4–10.7	4.7±1.3	5.1±1.1	5.7±1.0	6.4±1.1
	Crown base height variation (SD)	3.5±1.7	0.6–10.2	2.6±1.5	3.0±1.4	3.6±1.3	4.8±1.5
ASMRb	Mood disorders, women	118±33	40–221	119.3±27.6	121.9±26.0	121.0±33.8	110.1±40.0
	Mood disorders, men	122±34	40–241	129.0±29.8	126.8±26.9	124.3±32.4	109.2±42.6
	Cardiovascular disease, women	160±45	40–349	176.4±44.2	166.8±39.0	155.0±44.9	142.4±45.8
	Cardiovascular disease, men	151±40	40–334	152.7±33.5	156.2±33.3	152.7±41.8	142.4±47.7

Note: %LMIC, percentage foreign-born inhabitants from low- and middle-income countries; %Low education, percentage primary-educated-or-lower inhabitants; %Unemployed, percentage unemployed inhabitants; ASMR, age-standardized medication sales rate (prescriptions per 1,000 inhabitants); LiDAR, airborne light detection and ranging; Q, quartile; SD, standard deviation; SES, socioeconomic status.

The tree characteristics were derived from the 2012 winter airborne LiDAR data.

The ASMR data were collected between 2006 and 2014.

The accuracy of individual crown delineation was high: 91% of the delineated tree crown polygons contained a single tree stem, indicating a good agreement (Figure 3). In the 604 census tracts, 309,757 trees were included for the analysis. The average number of trees was 513 in total and 24/ha (IQR: 21.8) and the average crown volume was /ha (IQR: 1.4). Variation for tree height (expressed by SD) averaged and for crown base height . The geographical distribution of the number of trees and crown volume are presented in Figures S1 and S2, with the lowest numbers and the smallest volumes (lower 10 percentile values being 8 trees and /ha, respectively) found in the city center and along the canal and railway zone that bisects the city from northeast to southwest. The highest densities (85–142/ha) and crown volumes (/ha) were recorded in Sonien Forest along the southern edge of the region, but these census tracts are uninhabited and did not contribute to our models. The highest densities and crown volumes of inhabited census tracts were found along the perimeter of the study area (upper 10 percentile values being 50 trees and /ha, respectively).

Figure 3.

Subset of individual tree delineation results and LiDAR canopy height model and accuracy of individual tree segmentation. Note: LiDAR, airborne light detection and ranging.

Average medication sales for mood disorders were 118 sales/1,000 inhabitants for women and 122 for men. For cardiovascular medication, average sales were 160 /1,000 inhabitants for women and 151 for men. The geographical distribution of the medication sales for mood disorders and cardiovascular disease are presented in Figures S3–S6.

Associations between Urban Trees and Health

In unadjusted, single-factor models, both tree density and crown volume were associated with lower overall average medication sales. Higher tree density (an increase by one IQR of 21.8 trees/ha) was associated with 12.9% lower overall average medication sales [() and 95% confidence interval (CI): , )] and higher canopy volume (an increase by one ) with 16.3% lower overall average medication sales [ (95% CI: , )]. Higher tree densities were strongly correlated to higher crown volumes (Pearson , ). In a two-factor model for overall mean medication sales, higher crown volume was associated with lower medication sales [ (95% CI: , )] and higher tree density with higher medication sales [ (95% CI: 0.351, 0.408)]. These contrasting associations were conserved after adjusting for SES: higher crown volume (an increase by one ) was associated with 27.7% lower medication sales [ (95% CI: , )] and higher tree density (an increase by one ) with 26.3% higher medication sales [ (95% CI: 0.228, 0.298)].

The results of the detailed models for both medication types and sex groups separately are presented in Figure 4, and the parameter estimates for these associations are presented in Table S1 for mood disorder medication and in Table S2 for cardiovascular medication. Higher canopy volume (an increase by one ) was associated with 34% lower mood disorder medication sales [women, (95% CI: , ); men, (95% CI: , )] and with 21–25% lower cardiovascular medication sales [women, (95% CI: , ); men, (, )]. Conversely, higher stem density (an increase by one ) was associated with 28–32% higher mood disorder medication sales [women, (95% CI: 0.284, 0.361); men, (95% CI: 0.243, 0.319)] and with 20–24% higher cardiovascular medication sales [women, (95% CI: 0.169, 0.236); men, (0.206, 0.274)]. In addition, we found consistent but weak associations of crown structural variation (SD of crown base height) with medication sales for both medication and sex groups. For tree height variation, the associations were inverse but weak for all groups except for mood disorder medication sales among women, where no association was detected. Apart from minor variations in the magnitudes of the associations, no differences between men and women were found in the associations between tree exposure and medication sales.

Figure 4.

Associations between the number of trees, tree traits, and annual sales of medication commonly prescribed for cardiovascular disease and mood disorders in adults 19–64 years of age in Brussels, Belgium, at the level of census tracts (CTs, ). Parameter estimates and 95% confidence intervals (CIs) are from generalized log-linear models for age-standardized rates of medication sales in Brussels for the years 2006 to 2014. Models were adjusted for socioeconomic status at census tract level (percentage immigrants from low- and middle-income countries, percentage unemployed inhabitants, percentage inhabitants with only primary education).

In all stratified models, the associations of tree density and crown volume with medication sales were conserved (Tables S3–S6). In the analyses stratified by population density, the positive associations with stem density and the negative associations with crown volumes were stronger in less densely populated census tracts (), both for mood disorder and cardiovascular medication sales (Table S3). In the analyses stratified by unemployment and low education, associations were stronger in census tracts with low SES (, ) for mood disorder medication; for cardiovascular medication, associations were stronger in census tracts with high SES (, ) (Tables S4 and S5). In the analyses stratified by the proportion of inhabitants from LMIC, the estimates were comparable to those of the main models (Table S6).

Discussion

Main Findings

We investigated the association of sales of medication commonly prescribed for cardiovascular disease and mood disorders in Brussels, Belgium, with residential exposure to urban trees, adjusting for SES. Single-factor models using aggregated medication sales and exposure data showed that higher stem densities and higher crown volumes are both associated with lower medication sales. However, results from models that simultaneously included crown volume and tree density suggest that living in areas with high crown volumes and low stem densities may be more beneficial to adult cardiovascular health and mood disorders than living in areas with a higher amount of trees with smaller crowns.

Comparison with Other Studies

Previous studies using LULC data or street photography-derived tree canopy cover or inventory of street trees have found associations between urban trees and better cardiovascular and mental health (Akpinar et al. 2016; Astell-Burt and Feng 2019, 2020b; Donovan et al. 2013; Dzhambov et al. 2019; Egorov et al. 2020; Wang et al. 2020). Our findings that larger crown volumes were associated with lower mood disorder and cardiovascular medication sales are consistent with these earlier studies. Other studies, however, did not report such associations (Browning and Rigolon 2018; Jarvis et al. 2020). In addition, two earlier studies conducted in Brussels did not observe the expected benefits of tree/green space exposure on cardiovascular and mental health. Using the CORINE Land Cover data set and medication sales data, Aerts et al. (2020) observed an inverse association between reimbursed cardiovascular medication sales and tree/forest cover at the levels of census tracts in the entire country ( census tracts) but not in a separate model for Brussels ( census tracts), most likely because the resolution of the land cover data set was too coarse to adequately capture green space exposure in the heavily urbanized region of Brussels. In their study on the impacts of urban greenness on self-reported mental health ( participants) in Brussels, Pelgrims et al. (2021) quantified three more detailed indicators of urban greenness exposure: a) vegetation coverage derived from Google Street View panorama at the residence level; b) street tree density and vegetation cover ( on each side) from inventory data sets at the street level; and c) NDVI values at the neighborhood level. However, only street tree density was found to be inversely associated with sleeping disorders, whereas other urban green exposure indicators were not associated with any of the investigated mental health indicators. The individual trees in private areas, such as backyards and gardens, are most likely to be regularly experienced by urban residents (Cox et al. 2019) but are hard to represent in LULC maps in densely built zones where buildings still dominate the landscape as seen from 2D imagery (Becker et al. 2019). Moreover, the 2D tree cover derived from the LULC maps and street photography are less comprehensive for the actual exposure to trees than 3D indicators (Jiang et al. 2017). The lack of the associations in these two studies may thus be caused by the coarse representations of exposure to urban trees/green spaces, inadequately capturing numerous small green spaces and individual trees that are present in the study area and providing inaccurate information on tree dimensions. By identifying all trees in the study area and estimating crown volumes, our study improved the quantification of tree exposure and was able to observe health benefits of trees.

However, in two-factor models that included both crown volume and tree density, tree density was positively associated with medication sales. This is in contradiction with earlier observed associations between tree density per linear kilometer of street and antidepressant prescription rates reported by Taylor et al. (2015) and between street tree density within around the home and antidepressant prescriptions reported by Marselle et al. (2020). In addition, the indicators of structural variation in the present study show inconsistent associations with medication sales (Figure 4; Tables S1 and S2). Although the importance of tree diversity to human health is increasingly recognized, the findings on the associations between tree diversity and human health/well-being are still mixed (Aerts et al. 2018; Methorst et al. 2021; Sandifer et al. 2015). By including tree density, tree structural diversity indicators, and crown volume in the same model, our results reveal that there may be trade-offs between several tree traits with respect to their impact on human health. Our results suggest that larger crown volumes concentrated in fewer stems may be more beneficial for cardiovascular and mental health than similar crown volumes divided over more stems, despite the strong association between stem density and total crown volume. Large and well-developed crowns are one of the key attributes that distinguish large, old trees from young trees (Lindenmayer and Laurance 2017). The beneficial stress-reducing effects of large trees may be enhanced by psychological pathways toward restoration, including restoration of attention (Collins et al. 2020; Marselle et al. 2020; Taylor et al. 2015). In a study on the restorative effects of different forests, participants experienced significantly higher restorative effects when visiting old-growth and mature commercial forests than when visiting urban recreation and young commercial forests (Simkin et al. 2020). In urban environments, large trees may therefore help to reduce stress-related conditions, such as mood disorders and cardiovascular disease by, for instance, contributing more to indirect nature experiences than smaller trees (Cox et al. 2019) or by providing keystone habitat for birds (Le Roux et al. 2014; Pena et al. 2017; Stagoll et al. 2012), which in turn improve well-being via bird sightings and sounds (Buxton et al. 2021).

Strengths and Limitations

A major strength of this study is that we have retrieved 3D tree characteristics (most importantly, accurate tree crown volumes) for all the trees in the study area using very precise remote sensing–based data to determine accurate residential exposure to urban trees, which is a major improvement compared with previous research (e.g., Egorov et al. 2020; Marselle et al. 2020; Wang et al. 2020). The present study is also, to our knowledge, the first to specifically demonstrate health benefits of large trees in urban environments, indicated by the associations between crown volumes, stem densities, and medication sales data. This is of global conservation importance, because large, old trees in urban environments are expected to decline globally under current management practices (Le Roux et al. 2014). Nevertheless, our research may have a number of limitations.

The ecological study design is prone to ecological fallacy. Health effects of green space observed at the aggregated (census tract) level may not be observed at the individual level (Bixby et al. 2015; Mitchell et al. 2011; Richardson et al. 2012). Besides exposure to trees and SES (included as confounders in our models), other factors such as age, chronic conditions, and changes in prescription over time could also be strong predictors of medication sales (Webster 2007). The observed associations between exposure to trees and medication sales may be mediated by intergroup variations in age structure and incidence of chronic diseases, leading to inconsistent associations observed at the aggregated level and at the individual level. However, by using census tracts as the level of aggregation, we believe these intergroup variations have been minimized.

In addition, there may be a mismatch between the timing of exposure and medication use because of, for instance, seasonality. Seasonality of urban green has been recognized to affect participation in recreational physical activities (Villeneuve et al. 2018) and this may result in seasonal variations in associations between exposure to trees, human health, and therefore medication sales. We were not able to analyze the associations by season because information on the timing of medication use was not available in our data set and because our assessment of exposure to trees did not take into account potential seasonal variations in tree crown volumes. Future research might explore how the frequency and the timing of medication use would affect the tree–health associations.

The exposure measurements may contain uncertainties. Residential exposure to green space may differ from true exposure because it ignores dynamic time–activity patterns that include home time, school or working hours, recreation, and time in traffic and which may significantly impact health outcomes (Stas et al. 2021). The imperfect temporal match between the LiDAR data (2012) and medication sales data (2006–2014) may also lead to exposure measurement error owing to the changes in trees over time. One of the solutions to this issue is using multi-temporal LiDAR data, which is increasingly promising given that LiDAR data costs are decreasing and data processing is becoming more standardized (Zhao et al. 2018). In addition, like other tree identification and delineation algorithms, our algorithm was inevitably accompanied by issues, such as the misidentification of trees (e.g., building edges) and the over- and undersegmentation of tree crowns (Zhen et al. 2016). In addition, although the estimated crown volume using a 3D convex hull algorithm largely relies on the amount and distribution of branches that does not change much between summer and winter, for deciduous trees it may be underestimated owing to the use of winter (or leaf-off) LiDAR data (Ørka et al. 2010). We acknowledge that these issues may affect the tree characteristics aggregated at the census tract level. However, the building edge objects had a relatively small volume; the segmentation errors may slightly over- or underestimate the amount of trees (Figure 3) but did not bias the aggregated tree crown volume; and a recent study has shown that leaf-off LiDAR data were more successful in describing crown size diversity (Davison et al. 2020). Therefore, we believe these uncertainties did not meaningfully affect the tree–health associations observed in the present study. In addition, perfect identification and segmentation of individual trees in a highly heterogeneous environment with diverse tree species and local contexts (e.g., urban areas) is challenging and beyond the scope of this study. We acknowledge the fact that timing mismatch between medication sales data (2006–2014) and socioeconomic background data (2001) may have caused additional bias, but at the time of analysis this was the best possible match based on the availability of group-level socioeconomic and medication sales data sets.

Finally, it is important to note that lower medication sales may be associated with higher mortality. This may be the result of untreated cardiovascular disease and mood disorders, leading to higher mortality by cardiovascular disease and suicide. However, in a study on the associations between residential green space and mortality for residents in the five largest cities in Belgium (including Brussels), higher exposure to residential green space was associated with decreased nonaccidental mortality (Bauwelinck et al. 2021). Therefore, we argue that the lower medication sales in areas with higher exposure to green space in this study was most likely associated with the lower incidence of cardiovascular disease and mental disorders.

Conclusions

Based on aggregated health data and comprehensive 3D tree data, we found that both tree density and tree crown volume are inversely associated with medication sales for cardiovascular disease and mood disorders. However, results of models that evaluate several tree trait exposures simultaneously suggest that living in areas with large tree crown volumes divided over relatively fewer stems may be more beneficial for adult cardiovascular and mental health than living in areas with a similar crown volume divided over a higher number of trees with smaller crowns. Large tree crowns may reduce physical and mental stress more efficiently because the reduction of both heat and air pollution depend on leaf area, which is higher in large tree crowns. Psychological effects and indirect nature experiences provided by large trees, which are often old trees, may further strengthen the health impacts of these trees. Our results demonstrate that conserving large, old trees in urban environments supports not only the conservation of biodiversity but also human health.

Click here for additional data file.

Click here for additional data file.