Ray-Tracing-Based Numerical Assessment of the Spatiotemporal Duty Cycle of 5G Massive MIMO in an Outdoor Urban Environment

Featured Application: The presented numerical approach can be directly applied to the estimation of the compliance boundary of the antenna array base stations and the downlink human EMF exposure assessment in the networks served by such base stations. Abstract: In the near future, wireless coverage will be provided by the base stations equipped with dynamically-controlled massive phased antenna arrays that direct the transmission towards the user. This contribution describes a computational method to estimate realistic maximum power levels produced by such base stations, in terms of the time-averaged normalized antenna array gain. The Ray-Tracing method is used to simulate the electromagnetic ﬁeld (EMF) propagation in an urban outdoor macro-cell environment model. The model geometry entities are generated stochastically, which allowed generalization of the results through statistical analysis. Multiple modes of the base station operation are compared: from LTE multi-user codebook beamforming to the more advanced Maximum Ratio and Zero-Forcing precoding schemes foreseen to be implemented in the massive Multiple-Input Multiple-Output (MIMO) communication protocols. The inﬂuence of the antenna array size, from 4 up to 100 elements, in a square planar arrangement is studied. For a 64-element array, the 95th percentile of the maximum time-averaged array gain amounts to around 20% of the theoretical maximum, using the Maximum Ratio precoding with 5 simultaneously connected users, assuming a 10 s connection duration per user. Connection between the average array gain and actual EMF levels in the environment is drawn and its implications on the human exposure in the next generation networks are discussed. was implemented to determine the number of simultaneously active UEs. The results showed around 6 dB reduction for the 95th percentile of the maximum time-average BS gain as compared to the theoretical maximum, for high system utilization values.


Introduction
The much anticipated roll-out of the fifth-generation (5G) telecommunication networks brings about new challenges associated with limiting the exposure of the general population and workers to electromagnetic fields (EMF). One of the key universal features of the next-generation networks, shared among various 5G technologies, is the use of large antenna arrays at the base station (BS) side [1]. The radiation pattern of an antenna array depends on the amplitude and phase ratios of the array elements. By selecting the elements' amplitudes and phases in a specific way, a BS can produce directed "beams" in its far-field-the main lobes of the array radiation pattern. This technique is In [7], the approach of the authors of [6] was extended to a BS capable of MRT precoding (dubbed eigen-beamforming [7] or conjugate beamforming [8]). The 3GPP statistical model [9] was used with an urban environment, a single 64-element BS, and both indoor and outdoor UEs. Depending on the number of simultaneously served UEs, the 95th percentile of the compliance distance (which is proportional to the maximum time-averaged BS gain [5]) constituted 30% to 50% of the theoretical maximum.
This contribution builds upon and extends the approach proposed in [6,7] to more realistic precoding schemes and a more advanced Ray-Tracing (RT) channel model. The RT modeling yields spatially consistent channels, which depends on the environment geometry and UEs locations. For the MU-MIMO systems one important implication is a realistic UE DOD distribution, which governs the beam directions. The magnitude of the inter-UE channel correlation depending on the distance between the UEs is captured in the RT modeling-a factor that greatly impacts the BS pattern when using interference-canceling precoding schemes such as ZF. It has been shown that the RT method reproduces key parameters of measured Massive MIMO channels [10], whereas the state-of-the-art statistical model (WINNER-II) tends to underestimate the amount of correlation in the channels of closely spaced UEs [11]. Various other approaches to 5G channel modeling have been proposed in the literature and we direct an interested reader to the recent overview articles in [12,13].

Materials and Methods
In this section, we describe the environment model used in the RT simulations and methods for the results processing. An overview of the beamforming and precoding schemes is given.

Environment Model
The RT model consists of the environment geometry description and the transmitter-receiver (Tx-Rx) parameters. A geometrical entity is represented by the coordinates of its boundary faces (polygons) and the dielectric parameters (relative permittivity ε r and conductivity σ) of each face. A Tx (Rx) antenna is defined with its location, radiation pattern, and the carrier frequency. For every Tx antenna, rays are launched from its location in directions (nearly) uniformly distributed on a surface of a sphere centered at the Tx [14]. A ray is an abstraction which represents a flat wavefront, described with its (complex) EMF amplitudes and the propagation direction. The rays are propagated (traced) through the environment and their interactions (reflections, diffractions and transmissions), as well as path loss (PL) and time-of-flight, are tracked by the ray-tracer. If a ray passes sufficiently close to an Rx location, it is considered to be received and its state is recorded. The output of the simulation is a set of received rays for every defined Tx-Rx pair.
The RT output is site-specific, i.e., the channel between a fixed Tx-Rx pair depends on the surrounding geometry. To generalize the results a number of geometrical entities was generated stochastically, based on a few macroscopic parameters. Each realization of the environment geometry we call an environment sample [15]. Figure 1 presents one of the environment samples obtained from the model we describe below.
We simulate an outdoor urban macrocell bounded by a fixed flat square area 100 m by 100 m in size. Building blocks are represented by cuboids, width and length of which are sampled from a uniform random distribution in range from 15 m to 25 m. The height of a building block is drawn from a uniform random distribution in range from 5 m to 20 m. The buildings are positioned on a rectilinear grid, such that any two neighboring blocks are separated by exactly one empty grid cell. Rows of building blocks and straight lanes form a Manhattan-like urban city landscape. The spacing (lane) width is set equal to 10 m, 15 m, or 20 m randomly with equal probability. The dielectric properties of the cuboids model concrete material with ε r = 7, σ = 1.5 × 10 −2 S/m. The ground plane is assigned asphalt dielectric properties with ε r = 5.7, σ = 5 × 10 −4 S/m.
The locations and properties of the Tx and Rx antennas are fixed in the model. The simulations are performed at a single frequency of 3.5 GHz, foreseen to be heavily used in 5G networks.
The BS (Tx) is a rectangular 10-by-10 element array of vertically polarized half-wave dipole antennas with a half-wavelength uniform inter-element spacing (λ 85 mm at 3.5 GHz). The center of the Tx array is positioned at x = 1 m, y = 50 m and z = 25 m ( Figure 1). The BS height of 25 m correspond to the macrocell scenario in the 3GPP model [9]. In the following simulations, the BS coverage range spans from −60 • to +60 • in azimuthal (ϕ) and from 105 • to 135 • in polar angular directions (θ) in coordinate system shown in Figure 1. This is in accordance to the model used in [6]. To include most of the simulated ground-plane area within its coverage range, the BS array is tilted down by 30 • around the y-axis through its center.
A UE (Rx) is modeled as a single-terminal device equipped with a vertically-oriented vertically-polarized half-wave dipole antenna. The UEs are arranged on a regular rectilinear grid with 2 m spacing in x and y directions, at a height of 1.5 m above the ground-plane (z-axis). Only the grid nodes that fall within the cell and are located outside the building block interiors are kept in the simulation, which, on average, results in around 600 UEs per simulation. These UEs are used as potential active receivers in the analysis described in Section 2.3.
The density of the UE locations obtained in 25 environment samples is presented in Figure 2. The (ϕ, θ) coordinate system is shown in Figure 1. Figure 2 describes the UE locations in the DOD space, as viewed from the BS. The UE distribution averaged over the environment samples is symmetric with respect to the plane ϕ = 0, as a result of the matching (statistical) symmetry of the environment geometry. The UE density increases towards the upper cell edge (further away from the BS) as the polar angle of incidence decreases.
The RT simulation parameters were set as recommended in [16], limiting the environment interactions to up to 6 reflections, 1 diffraction, and 1 transmission.

MIMO Channel Matrix, Beamforming and Precoding
By s nk , n ∈ [1, 2, . . . , N], k ∈ [1, 2, . . . , K] we denote a set of rays between nth Tx and kth Rx points, with N and K being the total number of the Tx and Rx antennas, respectively. Then, the channel coefficientĥ kn is given by (see, e.g., in [17]) where u r is the voltage amplitude induced by the ray r at the Rx antenna terminal and τ r is the time-of-flight of the ray r. The MIMO channel matrixĤ is obtained by evaluating (1) for each Tx-Rx pair in the simulation. Columnsĥ k ofĤ are channel vectors to the Rx with index k. The distance from different UEs to the BS may vary significantly. Therefore, UEs will experience differences in PL in comparison to one another. Thus, channel equalization is performed by normalizingĤ column-wise [2] where || · || denotes the Frobenius norm.
Having the normalized channel matrix, the Massive MIMO precoding matrices W are given by [8] where (·) H denotes the Hermitian transpose and α is a real-valued normalization coefficient, chosen such that W has unit Frobenius norm. Furthermore, the codebook steering matrix W CB is constructed from the steering column-vectors b k ∈ C N as [18] In (4), d n is a distance vector from the BS array center to the nth element, c k is a unit vector in the codebook direction assigned to the kth UE, and (·, ·) denotes the dot product. Knowing the channel vectorĥ k , the beamforming direction c k is chosen as In (6), the maximization is carried out over the set of all beamforming directions c i supported by the BS. The beamforming direction vector c i is a unit vector in the direction (θ i , ϕ i ) of the ith beam center where θ and ϕ are the polar and azimuthal angles, respectively, in the spherical coordinate system depicted at Figure 1. The modeled system differentiates 32 beam directions in azimuth and 8 in elevation. The transmit vector t ∈ C N is obtained by multiplying the precoding or steering matrix by the vector of transmitted symbols s t = Ws.
As the EMF is further assessed in terms of the time-average root mean square (RMS) values, with no loss of generality we set all transmitted symbols to be real-valued. In addition, we assume that no per-user power management is implemented at the BS and equal share of transmit power is directed to each UE. Therefore, we define s = [ (9) with W given by (3) or (5) to satisfy the overall transmit power constrain.

Time-Average Antenna Array Patterns
An instantaneous array pattern is calculated as a sum of the patterns of its individual elements, weighted with the components of the transmit vector t. As all antennas in the BS array are identical dipoles, this gives where A dip is a half-wave dipole radiation pattern [19], and t n denotes the nth element of t. Here we do not account for the effect of mutual coupling in the antenna array, i.e., the modification of the free-space antenna element pattern by the currents in the neighboring elements.
In the far-field region of a BS, the incident EMF is proportional to the antenna gain in the direction where the measurement is preformed. As mentioned above, ICNIRP specifies [4] an EMF time-averaging interval T avg = 6 min for the human exposure assessment. At the same time, it is foreseen that in a typical scenario 5G DL traffic will be transmitted in short bursts (in the order of tens of seconds [6]), switching between sets of UEs that demand it at any given moment. If the served UEs are distributed uniformly enough within the cell, the BS would focus the transmission in many different directions over a sufficiently long time interval. Therefore, a realistic time-average BS antenna array gain is expected to differ significantly from the theoretical maximum one.
To quantify how much the time-averaged gain is reduced relative to the theoretical maximum we follow the approach proposed in [6,7]. We introduce a constant T-the duration of one connection ("drop duration" in [7] or "scheduling time" in [6]). We model a network in which independent sets of K UEs are served for time T in series with no overlaps. Then, the time-average BS array radiation patternÃ N,K m (θ, ϕ) is calculated as a weighted mean of the instantaneous patterns i produced during the averaging intervalÃ where m ∈ {CB, MRT, ZF} denotes the transmission precoding scheme used at the BS; N, K are the number of utilized antenna elements and the number of simultaneously served UEs, respectively, and ω i is the fraction of the averaging time during which pattern i was active, varying from 0 (not in the averaging interval) to T/T avg (fully inside the averaging interval). In the following, for convenience we choose T such that T avg is its integer multiple, then ω i = T/T avg . Next, the normalized gain G N,K m is given by where G N max is the maximum gain of an array of N elements. The maximum gain of a planar antenna array is calculated as a product of the maximum antenna element gain and the maximum array factor. The maximum array factor equals to the number of elements in the array [20]. Therefore, for an antenna array composed of identical half-wave dipoles, This value is further used in (12) as a normalization factor.
The RMS E-field strength at the location of the UE, to which G N,K m would be directed, can be estimated using a free-space approximation according to the following expression, where θ max is the polar angle of G N,K m , Z 377 Ohm is the impedance of free space, P is the BS total radiated power, and h BS = 25 m and h UE = 1.5 m are the BS and the UE height above the ground, respectively. For example, taking a BS with 64 antenna elements and 200 W nominal power [21], the highest achievable E RMS (G N,K m = 1) in the direction normal to the BS array plane (θ max = 120 • ) equals around 18.6 V/m, which falls below the ICNIRP reference level of 61 V/m [4].
To calculate G N,K m (T) for given m, T, N, and K, a numerical experiment is performed. Figure 3 shows a flowchart describing the procedure.
We studied configurations with 2-by-2, 4-by-4, 6-by-6, and 8-by-8 square sub-arrays selected from the center of the simulated 10-by-10 Tx array, as well as the complete array itself. These correspond to total array counts N of 4, 16, 36, 64, and 100 elements. Scenarios with K = 1, 2, 5, and 10 simultaneously active UEs were studied for each N. Connection duration T avg equal to 60 s, 10 s, and 1 s was considered for each N and K. In total, N env = 25 environment samples were simulated. Gain values obtained from the 25th sample were found to change the 95th percentiles of the time-averaged gain distributions by less than 1% (see Figure 5 and Section 4.2), which was accepted as a sufficient level of accuracy. Finally, in every environment sample N s = 100 time-averaged gain evaluations were performed, which amounts to 2500 evaluations of G N,K m (T) for each value of N, K, T and precoding scheme m. In the following section we present and discuss the distributions of G N,K m .
Select an environment sample.
Select K active UEs randomly out of all simulated Rx points. Calculate a vector of symbols s.
Select N antennas from the center of the BS array. Evaluate (1) and equalize (2) the channel matrix.
Repeat N avg = T avg /T times.
Repeat N s times with independently selected active UEs.
Repeat N env times with stochastically created environment samples.  Figure 4 shows the DODs of G N,5 CB (fist column), G N,5 MRT (second column), and G N,5 ZF (third column), for N = 4 (first row), N = 36 (second row), and N = 100 (third row), calculated with a connection time T = 60 s. The DOD of each time-averaged gain sample is depicted with a black circle in the (ϕ, θ) coordinate system. The circle size is proportional to G N,K m , and its opacity is proportional to the number of samples observed at the corresponding DOD. The background pseudocolor plot shows the sample-average ofÃ N,K m (θ, ϕ), illustrating the difference in the average beamwidths of the obtained patterns. G 4,5 ZF is undefined, as the condition N > K must be satisfied for HH H to be invertible, which is necessary to calculate G ZF according to (3). Its plot was therefore not included in Figure 4.  Figure 5 presents a CDF for K = 1 (black), K = 2 (red), K = 5 (blue), and K = 10 (green). In addition, the 95th precentile of each CDF is marked with a vertical dashed line of the same color. Table 1 lists the 95th percentiles for all possible combinations of the studied parameters. Additionally, the table cell background color saturation is proportional to its numerical value, ranging from white for zero to deep blue for one. This gives a visual cue to how different parameter combinations affect the normalized time-averaged gain. The circle opacity is proportional to the number of maxima found in the corresponding (ϕ, θ) direction. Left, center, and right columns show data for CB, MRT, and ZF transmission schemes, respectively. In the first, second, and third rows, scenarios with 2-by-2, 6-by-6, and 10-by-10 base station arrays are depicted, respectively. The ZF transmission with 2-by-2 BS array (N = 4) is undefined and was omitted. The dashed line depicts the cell boundary. The normalized time-and sample-averaged BS array patterns are shown in blue.

Array Patterns for CB, MRT and ZF
The left column of Figure 4 shows the DOD of the time-averaged gain observed when applying codebook precoding. The maximum value of the time-averaged array pattern was always found within the cell boundary. This was expected, as any instantaneous single-user codebook pattern has its main lobe pointing approximately towards an active UE, which was always situated within the cell. As a result of the linearity of (5), (9), and (10) with respect to the steering vector b, the 6 min time-averaged pattern is expressed as an average of the instantaneous patterns towards the UEs served during the averaging interval. The maximum of such averaged pattern was most likely to be found at the intersection of the instantaneous array patterns, i.e., somewhere within the cell. In the scenario with a 2-by-2 BS array (top-left), the maximum tends to be located around the cell center in azimuth. The reason for that is the low directivity of a typical pattern produced by an array of only 4 elements. For N = 36 (center-left), then the maxima distribution follows the average UE density peaks (Figure 2), with two clusters that correspond to the lanes between the building blocks, parallel to the x-axis (Figure 1). With 100 BS antenna elements (bottom-left), finally, the gain maxima closely follow the regions of high UE density (Figure 2), nearly covering the full azimuth range of 120 • .
The center column of Figure 4 shows the DODs of the time-averaged gain found using MRT precoding. Similarly to CB, when the antenna count is low (N = 4, top-center plot), the maxima tend to be concentrated around the cell center. As N increases, G MRT tends to be directed towards the regions densely occupied with the UEs with higher probability. However, MRT shows an increased spread of the gain DOD compared to CB. This can be attributed to the fact that unlike CB, which assigns a single beam per active UE, MRT superimposes a set of multiple beams with powers proportional to the contributions of the corresponding propagation paths to the total signal received by the UE. If a UE has a direct propagation path to the BS (i.e., LOS), the instantaneous MRT pattern is likely to have a global maximum in that direction (second strongest path-the ground reflection, if present, being orders of magnitude weaker). In case the target UE resides in a shadow region (NLOS), several propagation paths typically contribute comparable amounts to the total received signal, e.g., if a UE is obstructed by a building, the main propagation mechanisms that make the connection possible are over-the-rooftop diffraction and reflections from the walls of the surrounding buildings. As a result, the time-average pattern maximum is sometimes found outside the cell boundary, as can be seen on the plots showing G 36,5 MRT (center-center) and G 100,5 MRT (bottom-center) in Figure 4. This effect is even more pronounced when ZF precoding was used, although the underlying reason is different. ZF minimizes interference between the target UEs by canceling the transmission via shared paths. In scenarios with a large number of spatially correlated UEs, a portion of the total transmit power dedicated to fighting interference may exceed that of the intended signal, i.e., an instantaneous ZF pattern can have higher gain in its side-lobes than in the beams intended to reach the UEs. Such effect is observed in the DOD distribution of G 36,5 ZF (center-right) in Figure 4. In the areas with the highest UE density, where both CB and MRT generally produced their time-average gain maxima most often, ZF showed very few time-average array pattern peaks. The ZF precoding efficiency generally increases with the N/K ratio [2]. For N = 100 the gain distribution (bottom-right) was similar in shape to what was obtained using MRT, although the spread of the gain locations noticeably exceeded both MRT and CB.

Normalized Time-Averaged Gain
The CDFs in Figure 5 compare the G N,K m values for parameter configurations that correspond to those shown in Figure 4. Increasing the BS antenna count N decreases the normalized time-averaged gain for all studied schemes, with other parameters fixed. Two factors are contributing to this effect. First, the normalization coefficient in (12) is proportional to N, which counteracts the increase in the absolute array gain. Second, with larger N the BS is capable of producing narrower beams, which are less likely to interlap in the DOD region within the cell, reducing the maximum of the average taken according to (11).
Decreasing the connection time T also decreases the time averaged gain, that is, G N,K m (T = 60 s) > G N,K m (T = 10 s) > G N,K m (T = 1 s) for any fixed m, N or K. This observation is explained by the fact that the more independent UE sets are served in the averaging time-span T avg (or, equivalently, the less the T value is), the closer G N,K m approaches the normalized average of instantaneous array patterns A(θ, ϕ, t) over the cell. Conversely, in the limit of a single UE served with T = T avg , as follows from (10), (11), (12), the time-averaged gain is the instantaneous BS pattern maximum. In this case, the CB beamforming realizes the maximum theoretical gain G max for the codebook directions coinciding with the maxima of the BS antenna element's individual pattern. This can be seen by substituting t = b k into (10), assuming that the beam center DOD satisfies When only a single UE is connected at a time (K = 1, shown in black in Figure 5), CB, MRT, and ZF show very similar distributions of the normalized gain. In fact, as can be seen from (3) in the degenerate case of K = 1, the matrix inverse of HH H becomes a reciprocal of a squared channel coefficient magnitude. The ZF and MRT formulations are then equivalent, with appropriately chosen normalization coefficients α in (3). The minor discrepancy in Table 1 between the MRT and ZF (<1%) for K = 1 is due to the numerical round-off error propagation. The difference between G N,1 MRT and G N,1 CB , gradually increases with increasing N and decreasing T. Both G 4,1 MRT and G 4,1 CB are decreasing monotonously from around 0.85 for T = 60 s to around 0.66 for for T = 1 s (see Table 1). At N = 100, G N,1 m drops rapidly to around a half of that ( 0.42 for CB) for T = 60 s, around a third ( 0.22 for MRT) for T 10 s, and less than a quarter ( 0.15 for MRT) for T = 1 s.
Increasing the number of simultaneously served UEs K was found to decrease G N,K MRT and G N,K ZF for any fixed N and T. Increasing K decreased G N,K CB only for larger BS arrays (N ≥ 16) for T = 60 s and led to its increase with any N for shorter connection time values. For N = 4, the CB time-averaged gain closely approached the theoretical maximum for any T and K ≥ 5 (G N,K CB ≥ 0.9). This indicates that smaller BS antenna arrays implementing CB beamforming offer little to no benefit in terms of human EMF exposure reduction.
In a realistic usage scenario, we take T 10 s [6]. If the BS is equipped with 64 antenna elements, p 95 of G 64,K CB is just above 0.2 (around 7 dB reduction) for any K. This is in agreement with the results in [7] obtained with a similar configuration in the outdoor macrocell environment. Direct comparison with [6] is not possible, as in that case the UE count was varied during the averaging time. Increasing the averaging time T to 60 s increases the 95th percentile of G 64,5 m to 0.35 (4.6 dB) for CB, 0.30 (5.2 dB) for MRT, and 0.21 (6.8 dB) for ZF.
Adding more BS antenna elements while using either the CB or MRT scheme, does not decrease G N,K m significantly. Their lowest 95th percentile values were observed for T = 1 s: G 100,5 MRT G 100,10 MRT 0.14 and G 100,1 CB G 100,2 CB 0.13. The 95th percentiles of G N,K ZF were lower than the respective G N,K MRT and G N,K CB values for all K ≥ 2. This difference was larger with for larger K and shorter connection time T. In the realistic scenario with T = 10 s, the 95th percentiles of G 64,K ZF were equal to around 0.19 (7.2 dB), 0.13 (8.9 dB), and 0.1 (10 dB) for K = 2, 5, and 10, respectively. These values are nearly two times lower that the MRT and CB schemes demonstrated in the same parameter configurations. The lowest p 95 was found with G 100,10 ZF 0.06 (12.2 dB), which is around one third of the corresponding CB value, and more than two times lower than the minimum for the MRT or CB scheme.
As the time-averaged gain is directly related the average E RMS measured at some location in the cell. The far-field instantaneous E-field magnitude is proportional to the square root of the antenna gain. Therefore, the time-averaged E RMS is reduced at least in proportion to the square root of the time-averaged gain G N,K m , relative to the E RMS estimate based on the maximum achievable gain G N max . In a scenario with N = 64, K = 5, and T = 10 s, this leads to the E-field reduction in 95% of the observations by at least a factor of around 2.1 (3.2 dB), 2.2 (3.4 dB), and 2.8 (4.5 dB) for the CB, MRT, and ZF schemes, respectively, compared to the theoretical maximum. The theoretical maximum gain value was never reached in samples with K ≥ 5.

Conclusions
We presented a numerical approach that utilizes the RT method to model a time-averaged array gain of a 5G BS operating in a macrocell outdoor urban environment. The RT approach provides a more realistic signal propagation and user spatial correlation properties compared to analytical and stochastic approaches. In a realistic scenario, with a BS consisting of 64 antenna elements that serves 5 UEs simultaneously and a 10 s per user connection duration, 95% of the 6-min time-average gain observations fell below 0.22 (more than 6.6 dB reduction), 0.20 (7 dB) and 0.13 (8.9 dB) of the theoretical maximum, using codebook, maximum ratio transmission, and zero-forcing schemes, respectively. With user connection duration of 60 s, the corresponding 95th percentiles increase to 0.35 (4.6 dB), 0.30 (5.2 dB), and 0.21 (6.8 dB), respectively. In all studied scenarios, increasing the BS element count decreased the normalized time-average gain. With the MRT and ZF transmission schemes, lower time-averaged gain was always observed when the number of multiplexed UEs was increased. With the CB beamforming that was the case only for larger BS arrays. In all multi-user scenarios, the ZF yielded the lowest p 95 values of the normalized time-average gain (0.06 or 12.2 dB reduction with 100 BS antennas and 10 UEs), which is more than two times lower than any other studied precoding scheme. Funding: Piet Demeester thanks the ERC for his advanced grant 695495 "ATTO: A new concept for ultra-high capacity wireless networks". Arno Thielens is an FWO senior postdoctoral fellow. Sam Aerts is an FWO postdoctoral fellow.

Conflicts of Interest:
The authors declare no conflict of interest.