Monitoring the setting of calcium sulphate bone-graft substitute using ultrasonic backscattering

TL;DR

This paper presents an ultrasonic backscattering method to monitor the setting process of calcium sulphate bone-graft substitutes in real-time, enabling differentiation between liquid and solid states during setting.

Contribution

The study introduces a novel ultrasonic backscattering technique for real-time monitoring of calcium sulphate setting, linking acoustic properties to material state.

Findings

Backscattering analysis can distinguish liquid and solid states.

Correlation of SOS and AIB identifies setting transition.

Threshold of 15 dB effectively estimates setting state.

Abstract

We report a method to monitor the setting process of bone-graft substitutes (calcium sulphate) using ultrasonic backscattering techniques. Analyzing the backscattered fields using a pulse-echo technique, we show that it is possible to dynamically describe the acoustic properties of the material which are linked to its setting state. Several experiments were performed to control the setting process of calcium sulphate using a 3.5-MHz transducer. The variation of the apparent integrated backscatter (AIB) with time during the setting process is analyzed and compared with measurements of the speed of sound (SOS) and temperature of the sample. The correlation of SOS and AIB allows to clearly identify two different states of the samples, liquid and solid, in addition to the transition period. Results show that using backscattering analysis the setting state of the material can be estimated with a threshold of 15 dB. This ultrasonic technique is indeed a first step to develop real-time monitoring systems for time-varying complex media as those present in bone regeneration for dental implantology applications.

Figures

Tables

Table 1. Table 1: Fitted parameters for the phenomenological model given by Eq. ( 2 ).

	Transition asymmetry	Transition speed	Transition time	Initial state SOS	Final state SOS
	$\beta$	$\gamma$ (min-1)	$t_m$ (min)	$c_1$ (m/s)	$c_2$ (m/s)
Sample 1 (E1)	2,1	0,346	20,7	1477,2	2731,5
Sample 2 (E2)	1,9	0,341	20,8	1399.7	2490.9
Sample 3 (E3)	1,4	0,327	19,9	1553,9	2919,5
Sample 4 (E4)	1,7	0,367	20,1	1627,2	3049,8
Mean	$1,8\pm 0.3$	$0,345\pm 0,016$	$20,4\pm 0.4$	1514.5$\pm$98.0	$2858.1\pm 155.5$

Equations

$$c=\frac{2L}{\Delta t},$$

$$c(t)=c_1+(c_2-c_1){\tan }^{-1}{\left[\frac{\gamma }{\pi }(t-t_m)+\frac{1}{2}\right]}^{\beta },$$

$$ABTF=10{\log }_{10}\left(\frac{P_s(f)}{P_{\mathrm{ref}}(f)}\right),$$

$$AIB=\frac{1}{\Delta f}\underset{f_1}{\overset{f_2}{\int }}ABTF(f)df,$$

Full text

I Introduction

Dental implants require sufficient bone volume to achieve stability and durability. The loss of tooth affects the quality and thickness of the bone Bodic et al. (2005). In cases where bone volume is insufficient, guided bone regeneration (GBR) techniques are appliedBodic et al. (2005). In particular, GBR by maxillary sinus floor augmentation procedures are used since 50 years, and were first proposed by Tatum in 1974 Tatum (1986). In this operation, a window in the lateral area of the gingival is performed, the Schneider membrane is lifted, and a bone graft inserted. The window is coated with a pericardium membrane. Full bone regeneration process has a typical healing duration of approximately 6 months. The procedure sometimes failAnnibali et al. (2008) and different complications may appear postoperatively, e.g. soft tissue growth or failed osseointegration Acero (2011). Odontologists make use of X-ray radiological techniques in order to evaluate the bone regeneration process Memon et al. (2010). However, monitoring the full regeneration process using these techniques imply potential risks, such as excessive exposure to ionizing radiation Memon et al. (2010); Ludlow, Davies-Ludlow, and White (2008). In this sense, non-invasive ultrasonic techniques are desirable to monitor bone-regeneration during the healing time since they are non-ionizing.

Ultrasound has been used fruitfully to characterize the physical properties of bone. Technology for bone density estimation is nowadays mature Abendschein and Hyatt (1970); Rossman et al. (1989). Many commercial devices exist in the market used for the diagnosis of osteoporosis, e.g., for calcaneus or phalanx Langton, Palmer, and Porter (1984); Langton and Njeh (2008). Acoustic and elastic parameters are correlated with the physical properties of bone for diagnostic purposes. Among them, the speed of sound (SOS) is the most spread (Laugier and Haïat, 2011, Chapter 3). It consist in estimating the effective travelling velocity of a pulsed excitation, associated to the group speed of longitudinal waves in a given frequency band. Bone tissue is a heterogeneous biphasic medium composed of a multiple-scale porous viscoelastic matrix, soft-solids and viscous fluids Laugier and Haïat (2011). Acoustic propagation through these complex media is therefore strongly dispersive and SOS, given as a single value, should be interpreted as an effective parameter rather than an intrinsic physical parameter locally related to the mechanical properties of the bone tissue. A simple way to measure the SOS of bone is to estimate the time-of-flight difference between two received echoes using a sample of known thickness.

In addition, ultrasonic waves travelling in bones are strongly scattered in a diffuse way due to their complex heterogeneous micro-structure. In this way, the interference of multiple scattered waves from internal micro-structures can be received at source position. During last two decades, different techniques based on the backscatter energy evaluation have been proposed for bone characterization Wear and Garra (1998); Chaffaı et al. (2002). Ultrasonic backscatter energy can be linked to the physical properties of bone tissue such as bone volume fraction or bone mineral density Laugier and Haïat (2011); Padilla et al. (2008); Karjalainen et al. (2009), enabling quantitative ultrasound techniques for bone characterization.

Many ultrasonic backscattering techniques are based on measurements of the apparent backscatter transfer function (ABTF)Laugier and Haïat (2011). It represents the backscattered power from the sample corrected for the frequency response of the measurement system. Two parameters without explicit dependence of frequency are obtained from ABTF. The apparent integrated backscatter (AIB) is one of the most spread. It consist of a measure of the frequency-averaged (integrated) backscattered power contained in some portion of a backscattered ultrasonic signal Hoffmeister et al. (2006). The frequency slope of apparent backscatter (FSAB) is determined by fitting a line in the frequency domain to the spatially averaged backscattered ultrasonic signals over the analysis bandwidth Hoffmeister, Holt, and Kaste (2011). Both AIB and FSAB are determined in a logarithmic scale and are referenced to the measured reflected energy of a reference flat panel, e.g., a steel sample. In order to remove this reference, Hoffmeister proposed in Hoffmeister et al. (2012) a set parameters that compare the backscatter energy from two different parts of the sample. The mean of the backscatter difference spectrum (MBD) is obtained by frequency averaging the difference between the backscatter spectrum over the analysis bandwidth from two different windows of the ultrasonic signal.

During bone regeneration process, its mechanical properties are modified from the initial bone graft to the consolidated bone. A change of state is produced as a consequence of osseointegration. Thus, a solution for detecting possible anomalies during the healing of bone in GBR applications consist of monitoring the ultrasonic properties of tissue during time. In this sense, ultrasound techniques have been applied for characterization of slowly time varying materials. One relevant example are cements, that have been characterized by ultrasonic methods for more than 80 years, see e.g, Naik, Malhotra, and Popovics (2003); Philippidis and Aggelis (2005). In particular, pulse-echo ultrasonic techniques can be used to characterize the time dependency of SOS and the acoustic impedance during their setting processKeating, Hannant, and Hibbert (1989); Sayers and Grenfell (1993); Boumiz, Vernet, and Tenoudji (1996). Cement starts to set when mixed with water and a series of hydration chemical reactions are produced. The constituents slowly hydrate and the mineral solidifies. In particular, the typical setting duration time for bone-graft materials based in calcium sulphate ranges from 30 minutes to 60 minutesCarlson et al. (2003); Vlad et al. (2012). Same as cement setting, the GBR is a process in which the mechanical (and acoustical) properties of the material vary with time and a change of state is produced. The advantage of using a bone graft substitute as a phantom for monitoring GBR is that the duration of the process is reduced from 6 months-time of healing, to less than an hour in a controlled environment.

In this work, we propose a proof of concept to monitor the GBR by sinus augmentation using ultrasonic backscattering methods. Thus, the time-varying acoustical properties of the material under inspection are obtained using an echo-impulse technique. In particular, we use a phantom of calcium sulphateFinkemeier (2002); Wang and Yeung (2017), which has already been used for GBR in the past De Leonardis and Pecora (2000); Thomas and Puleo (2009). The change of state from a viscous fluid to a porous solid during the setting mimics the GBR in a much shorter duration that can be studied under laboratory conditions. We report that ultrasonic backscatter energy can be linked to the physical properties of tissue. In GBR applications the transmission mode using a pair of ultrasonic transducers is not allowed as a receiver transducer can not be placed behind the area under inspection. Thus, we propose the use of a single sensor improving the accessibility of medical instrumentation for tissue characterization.

II Methods

II.1 Phantom preparation

Calcium sulphate type IV (Ventura Pinkmod) is the synthetic bone substitute used to perform the experiments. In particular, in its powder state the material is calcium sulfate hemihydrate (\chCaSO4 . 1/2 H2O) in its $\alpha$ form. When the hemihydrate is mixed with water, calcium sulphate dihydrate is formed in a moderate exothermic reaction asThomas and Puleo (2009):

\ch

CaSO4 . 1/2 H2O + 1_1/2 H2O -¿ CaSO4 . 2 H2O + Heat.

In this way, during the setting process the material shows a change of phase from a viscous liquid to porous solid. The total setting time is around 20 minutes. To study repeatability, four experiments were performed with samples composed of in 23 ml of water and 50 g of calcium sulphate. The powder was dissolved in water at 20 $\pm 0.2^{\circ}$ C and mixed up to get homogeneous samples. The material was deposited in a cylindrical plastic container with a diameter of 36 mm to provide samples of 18 mm of thickness.

II.2 Ultrasonic measurement

An Olympus V382 transducer with a nominal central frequency of 3.5 MHz and a bandwidth of 2.34 MHz (-6 dB) was used as emitter and receiver using a pulse-echo technique. Figure 1 shows the block diagram and a photograph of the experimental setup. The transducer was placed in contact to the bottom of the recipient with the sample. Vaseline was applied between the transducer and the container to enhance acoustic coupling, and between the container and the sample to avoid decoupling during the solidification process. An ultrasonic pulser-receiver (Panametrics Model 5072) was used for emitting and receiving the acoustic pulses. The received signal was digitized with a data acquisition platform (Red Pitaya) at a sampling frequency of 125 MHz. For each experiment, a total of 1600 acquisitions were performed, each one every 2 seconds. The total duration of each experiment took 53 minutes, enough to cover the full setting process of the samples. Temperature was registered with thermistor probe (Tinytag Temperature Logger TK 4023) located externally in contact with the container. Measurements were post-processed to provide the acoustic parameters during the setting process.

II.3 Speed of sound

The speed of sound (SOS) was calculated as the distance ratio between two times the thickness of the sample and the time of flight between the first two echoes from the end of the sample as

[TABLE]

where $L$ is the thickness of the sample and $\Delta t$ is the time of flight, estimated as the times corresponding to the peak value of first two echoes. SOS was estimated during setting time, and a purely phenomenological model was fitted. The proposed model describes a smooth and asymmetric transition between two phases, fluid and solid, and is given by:

[TABLE]

where $c_{1}$ is the speed of sound in liquid state, $c_{2}$ is the speed of sound in solid phase, $\gamma$ models the speed of the transition from liquid to solid phase, $t$ is the setting time, $t_{m}$ is the half-time of the setting process and $\beta$ is related to the symmetry of the function. Using this model we can compare the setting process of several experiments.

II.4 Backscatter parameters

Two quantitative parameters based on the backscatter energy were considered in this work: the apparent backscatter transfer function (ABTF) and the apparent integrated backscatter energy (AIB). Both have been proven to be reliable and robust to characterize bone tissueHoffmeister et al. (2015). On the one hand, the ABTF is given by the transfer function between the backscatter signal and a reference signal as

[TABLE]

where $P_{s}(f)$ is the frequency-dependent power of a backscatter signal $p_{s}(t)$ using a time window of duration $\tau_{w}$, and $P_{\mathrm{ref}}(f)$ is the frequency-dependent power of the first echo from a reference reflector, e.g., steel plate Hoffmeister et al. (2015). On the other hand, the AIB is obtained as the frequency-averaged ABTF as

[TABLE]

where $\Delta f=f_{2}-f_{1}$, and $f_{1}$ and $f_{2}$ are given by the -6 dB bandwidth of the transducer. The reference pulse and its spectrum are shown in Figs. 2 (b,c), respectively. The corresponding values for the effective limiting frequencies are $f_{1}=3.4$ MHz and $f_{2}=4.6$ MHz.

III Results

III.1 Description of backscatter signals

The experimental signals of one experiment are shown in Fig. 2. First, in Fig. 2 (a) the A-lines are presented as a function of the setting time. Here, the amplitude of the ultrasonic signal registered by the transducer is represented in logarithmic scale normalized to the maximum value for visualization. The vertical axis represents the time of the experiment during the setting process (in minutes) and the horizontal axis represents the time of the ultrasonic signal (in $\mu$s). As both scales are very different ($\sim 10^{6}$), they are decoupled and each ultrasonic capture is assumed to be instantaneous.

From $t=0$ $\mu$s to $t=2$ $\mu$s, all ultrasonic signals have a maximum value of 0 dB corresponding to the emission pulse. Then, vertical grey traces appear corresponding to the received backscattered waves. First, when the experiment is started, the sample presents a viscous fluid state. At the beginning of the experiment time, at the top of Fig. 2 (a), the echo from the end of the sample is visible at $t=23$ $\mu$s. Up to about 15 min, the backscatter signals are stable and only small changes are visible, as shown by the straight vertical traces in the map. However, from about 18 min to 27 min a sudden and drastic change in the signal is observed: the line corresponding to the echo from the back of the sample is shift towards the beginning of the signal. This temporal compression of the signal is related to the increase of SOS in the material: the echo from the back is arriving sooner, as we will see after. Then, at 27 min the signals become stable again and, therefore, SOS should present a stable value at the end of the experiment. At 50 min the sample is physically a porous solid material.

Dashed lines in Fig.  2 (a) show the limits for the temporal window used for the post-processing of the backscattering. For their estimation, both the initial and final temporal values were corrected using the SOS in order to capture a constant spatial window. Figure 2 (f) shows in detail the tracking of lines during the transition from liquid state to solid state in the setting process. Here, the traces of the captured backscattered waves match the analysis window marked by the dashed white lines.

Backscatter parameters have been referenced to a signal obtained from a flat steel plate, shown in Fig. 2 (b). An example of A-scan signal registered by the transducer using the bone graft substitute sample (at $t=25$ min) is shown in Fig. 2 (d). Both the initial electrical contribution and the first echo caused at the end of the sample can be observed at $t=0$ $\mu$s and $t=15$ $\mu$s, respectively. In comparison, backscattered energy is very low and spreads all along the signal. We consider a window of duration of $\tau_{w}$ = 2.8 $\mu$s, i.e., 10 times the period corresponding of the central frequency of the transducer at the beginning of the experiment, according to Hoffmeister et al. (2006). The window was set initially centred at 8 $\mu$s. Note that this value is shifted at 25 mins (Fig. 2 (d)) to 5 $\mu$s due to the change of SOS with the setting time (see white dashed lines in Fig. 2 (a)). The spectra of the reference signal and the example backscattered signal are shown in Figs. 2 (c,e), respectively. Then, for each temporal signal ABTF is calculated using Eq. (3), as shown in Fig. 2 (f). Finally, AIB is estimated using Eq. (4) by averaging ABTF over the interest bandwidth, shown in Fig. 2 (f) as the shaded region.

III.2 Evolution of speed of sound

The evolution of the SOS with the setting time is shown in Fig. 3 (a). In addition, the evolution of the temperature and its variation, i.e., $\partial T/\partial t$, are shown in Fig. 3 (b). The initial value measured for the SOS corresponds to $c_{1}=1477.2\pm 50$ m/s close to the value of liquid water at same temperature (1495 m/s). Up to $t_{1}=17.32$ min, temperature is approximately stable and its rate of variation is small, about $0.1^{\circ}$C/min. However, after this time ($t_{1}$), the chemical reaction between the calcium sulfate and water suddenly increases. When temperature starts to increase small variations of SOS were already detected.

The maximum temperature variation is achieved at 21.6 min, corresponding to $1.9^{\circ}$C/min. This point corresponds to the maximum exothermic chemical activationLiu et al. (2003). Note this time almost corresponds to the time $t_{m}=20.7$ min given by the SOS model fit. After that moment, there is more consolidated material than calcium sulphate hemihydrate and water available to maintain the chemical reaction at this rate and, therefore, the exothermic activity is smoothly reduced. Then, at $t_{2}=27.6$ min the temperature reaches its maximum value, $35.5^{\circ}$C. The sample starts to cool smoothly due to heat diffusion processes and the lack of chemical activity. Once the change of state is completed, the sample is transformed from a viscous liquid to porous solid material. When the material is solidified, the SOS reaches a stable value close to $c_{2}=2731.5$ m/s. The behaviour of this synthetic bone substitute agrees the one reported for other cements measured using a emission/reception modeKeating, Hannant, and Hibbert (1989); Sayers and Grenfell (1993); Boumiz, Vernet, and Tenoudji (1996).

As expected, it can be observed that the change of the state of the material corresponds to an abrupt increase of the SOS. This particular smooth transition was measured in all samples, and can be described by three different phases: the liquid state, the transition and the solid state. The experimental SOS data was fitted to the phenomenological model given by Eq. (2). The parameters of the fitted model describe quantitatively the setting process. These parameters are listed in Table 1 for the 4 samples, showing similar values between experiments.

First, the parameter $c_{1}$ is linked to the physical SOS at the initial viscous-liquid state of the sample. The mean value of $c_{1}$ over the four experiments is 1514.5 $\pm$ 98.0 m/s. As the material at the initial process is mainly water with undissolved powder, the SOS measured for the sample is close to the SOS of the water. The final value of the SOS in the experiment is related to $c_{2}$, i.e, the SOS of the material in solid state. The mean value of $c_{2}$ is 2858.1 $\pm$ 155.5 m/s.

The transition time for the setting process is given by the parameter $t_{m}$. This parameter takes a mean value of $t_{m}=20.4$ min, very close to the peak of maximum exothermic chemical activity (21.6 min for sample 1). The factor $\gamma$ is related to the duration of the transition from liquid to solid state, its mean value is 1.8 $\pm$ 0.3 min*-1*. Finally, the factor $\beta$ is related to the symmetry of the function and takes a mean value of 0.345 $\pm$ 0.016. The dispersion of the parameters is low, and all experiments present similar transition curves. However, some discrepancies are observed in SOS values, that can be associated to uncertainness on estimating the thickness of the sample $L$.

III.3 Evolution of apparent integrated backscatter energy

In addition to SOS, energetic parameters also can provide information about the change of phase of calcium sulphate. In particular, we show in Fig. 3 (c) the AIB as a function of the setting time. Moreover, the amplitude of the 1st echo is also shown by the red curve in Fig. 3 (c). When the sample is in a viscous-liquid state, from $t=0$ to about $t_{1}$, the mean value of AIB is low and constant, taking a value of about $-15$ dB. However, the AIB estimations present high dispersion due to the low amplitude of the backscatter signals. Note that, in addition, the amplitude of the 1st echo is also very low. The weak scattering and the weak echo reflection are caused by the strong absorption and high homogeneity of the sample in its initial viscous-liquid phase. When the exothermic chemical reaction starts, both AIB and echo signal initially decrease. AIB takes a value down to -20 dB at 20 min.

However, as the reaction continues (after minute 17) both the scattering energy and the peak amplitude of the echo increase: the intrinsic absorption of the consolidated synthetic bone substitute is lower compared to its viscous-fluid phase. Thus, AIB progressively increases and its dispersion is reduced. During this time the transition of state in the setting process is produced; the mechanical and acoustic features of the sample are drastically modified and the microstructure of the porous solid material starts to consolidate, leading to a increase of the amplitude of the backscatter signal. The consolidated material present less absorption, and at the middle of the transition the echo from the back arrives with remarkable amplitude. However, as the porous solid consolidates its inner microstructure the backscatter energy increases. As a consequence, the amplitude of the 1st echo of the back is again reduced as waves reach the end of the sample with lower amplitude. After some small oscillations during the transition the AIB ends with a value of about -5 dB at 25 min.

Finally after 30 min, the sample is transformed to solid state and the AIB is smoothly increased from -5 dB to a plateau reaching a maximum value 0 dB at 50 min. Therefore, while the dispersion of the AIB is low and it reaches a stable value, the amount of the energy of the ultrasonic pulse that reaches the end of the sample is reduced due to scattering processes. This is caused by the final properties of the bone graft substitute that present high micro-heterogeneity and porosity. After this time, the setting process is completed and the material is mechanically and acoustically stable.

III.4 Correlation between of SOS and AIB

While at the initial and the end of the setting time AIB takes different values, AIB presents some oscillations during the setting period. These can lead to uncertainness if used to monitor the setting of the material. A deeper analysis can be performed by showing the correlation between AIB and SOS. Figures 3 (a-d) show this correlation for the four experiments. For each experiment, the insets show the corresponding histograms for SOS and AIB data, where two peaks corresponding to the initial and final states can be observed. While there exist a positive correlation between AIB and SOS, the data show some oscillations.

A concentration of experimental points can be observed in two zones, marked by arrows. On the one hand, when the material is in its viscous liquid phase is when the values of SOS are minimum. On the other hand, when the material is a porous solid both SOS and AIB present maximum values. As experimental data concentrates in these two limits, we can apply clustering techniques to identify them. A density-based spatial clustering of applications with noise (DBSCAN) techniqueEster et al. (1996) was employed here to discern between the data corresponding to the two concentration points and the transition. This give us the corresponding clusters marked in red and blue in Figs. 3 (a-d), helping to monitor the setting process: when a clusters is clearly populated, the sample has reached a stable phase. Note that using only the histogram for the AIB some of the signals corresponding to the transition point can present the same AIB than in the final state due to the oscillations of this parameter, see e.g. Fig. 4 (c).

IV Conclusions

The change of state of a bone graft substitute (calcium sulphate) has been monitored and characterized by ultrasonic backscattering analysis. Temperature, speed of sound (SOS) and apparent integrated backscatter energy (AIB) were measured during the setting process of the sample using a pulse-echo technique. The evolution of SOS shows a smooth transition between two limiting values, corresponding to the viscous-fluid and porous solid state of the calcium sulphate samples. An excellent agreement is found between the evolution of the SOS and a simple phenomenological model to describe the setting process. The experiments also show good repeatability.

In addition, the evolution of AIB is shown to be correlated with SOS. The apparent integrated backscatter energy is proven to be an appropriate parameter to describe the transition from a viscous-liquid state, with values around -15 dB, to a complex porous solid, with values around 0 dB. This relation confirms that the backscattering analysis with ultrasound is a powerful tool to monitor complex materials with time-varying mechanical properties using a single transducer and pulse-echo techniques.

This study represents a contribution to the use of ultrasound for the monitoring of GBR, which is necessary for the successful placement of dental implants. A synthetic bone substitute has been used as a bone phantom because its mechanical properties change during the setting process from a viscous liquid to a rigid porous material, as, under reasonable simplifications, occurs in GBR applications. However, many practical questions remains open and should be explored. These include the effect of the varying coupling of the transducer and the tissue between acquisitions, the mechanical/acoustical behaviour of real tissue during bone regeneration, the scattering properties of these unconsolidated and consolidated tissues, patient variability, or the impact of the intrinsic properties of real bones during regeneration (porosity, tortuosity, viscosity, …) in the estimated parameters. These issues should be explored in future studies.

Finally, it is worth to mention that one advantage of backscattering analysis with respect to speed of sound is that using backscattering techniques there is not needed to know the thickness of the sample. Note using a single transducer the estimation of the thickness of the sample is critical to accurately estimate SOS and, therefore, scattering techniques might be preferred. This is a significant benefit to envisaged implant applications in-vivo where the GBR can be monitored.

Acknowledgment

This work has been developed within the framework of the IVIO-UPV Chair, of the research project of the Universitat Politècnica de València “Desarrollo de tecnologías aplicadas al campo de la odontología” and is the result of the collaboration of researchers from the Universitat Politècnica de València and the Instituto Valenciano de Investigaciones Odontológicas (IVIO) in the field of dental implantology with ultrasonic technology. This work has also been supported by UPV and FISABIO through the project OSEODENT (POLISABIO 2018). The authors acknowledge financial support from European Union and Generalitat Valenciana through the European Regional Development Fund program (IDIFEDER/2018/022) and Agència Valenciana de la Innovació through Unitat Científica d’Innovació Empresarial (INNCON00/18/9). N.J. acknowledges financial support from Generalitat Valenciana through grant APOSTD/2017/042.