**A Versatile Setup for Time-Resolved Functional Near Infrared Spectroscopy Based on Fast-Gated Single-Photon Avalanche Diode and on Four-Wave Mixing Laser**

#### **Laura Di Sieno 1, Alberto Dalla Mora 1, Alessandro Torricelli 1,2, Lorenzo Spinelli 2, Rebecca Re 1,2, Antonio Pi**ff**eri 1,2 and Davide Contini 1,\***

Received: 15 April 2019; Accepted: 4 June 2019; Published: 10 June 2019

**Abstract:** In this paper, a time-domain fast gated near-infrared spectroscopy system is presented. The system is composed of a fiber-based laser providing two pulsed sources and two fast gated detectors. The system is characterized on phantoms and was tested in vivo, showing how the gating approach can improve the contrast and contrast-to-noise-ratio for detection of absorption perturbation inside a diffusive medium, regardless of source-detector separation.

**Keywords:** time-domain NIRS; null source-detector separation; brain; diffuse optics

#### **1. Introduction**

Among the fascinating applications of lasers, the non-invasive/non-destructive analysis of diffusive samples [1] (such as biological tissues, fruit, chemicals and wood [2–4]) has recently gained a lot of interest for the promising exploitations in biomedicine, in the food sector and in industrial fields. Photons injected in a diffusive medium undergo many scattering events due to refractive index mismatches at the micro and meso scale. Moreover, if light absorption is weak, like what happens in the so-called "therapeutic window" (600 to 1100 nm), where biological tissues are almost transparent to red and near infrared light, diffusing photons can easily be reemitted at the medium boundaries after having travelled deep in the tissue (e.g., >1 cm below skin in biological tissue in a reflectance configuration). The information carried by diffusing photons relates to the structural properties (linked to the scattering coefficient) and to the chemical properties (linked to the absorption coefficient) of the medium. For biomedical applications, the main chromophores of interest in the range between 600–850 nm are oxygenated and deoxygenated hemoglobin (O2Hb and HHb, respectively) [5]. For this reason, diffuse optical spectroscopy is a promising tool to assess brain and muscle oxygenation. In particular, functional near infrared spectroscopy (fNIRS) exploits the difference in absorption spectrum between HHb and O2Hb to recover their concentration over time, thus enabling the ability to see variation in haemodynamics, which can be a fingerprint of a brain activation thanks to the neurovascular coupling mechanism.

While the use of steady state light sources and detectors enables the so-called continuous wave (CW) approach to fNIRS, the time-resolved (TR) approach relies on the injection in the sample of short laser pulses (few tens of ps duration) and on the recording of the arrival time of re-emitted photons (distribution of time of flight, DTOF) at a given distance, the so-called source-detector distance (SDD). It has been demonstrated that the TR-technique provides a higher information content than CW, in particular when a single SDD is used [6]. Indeed, it is possible to recover the optical properties (absorption coefficient -μa- and reduced scattering coefficient -μs'), after fitting of the DTOF with a proper theoretical model [7]. Moreover, when using the TR technique in the so-called "reflectance geometry" (i.e., injection and collection put on the same side of the sample), the depth reached by photons is intrinsically encoded in their arrival time—the later the photon is re-emitted, the longer is its pathlength, thus increasing significantly the probability of the photon to have probed a deeper region [8]. It has also been demonstrated that the mean penetration depth of photons is independent of the SDD. For this reason, the null-SDD approach was proposed some years ago [9]. The use of short (or even null) SDD improves the confinement of photons, thus resulting in a higher spatial resolution and number of photons recorded, as well as a larger contrast at all times. To apply this approach, however, a fast-gated detector (i.e., a detector that can be turned on/off in a few hundred ps) is needed to prevent the detector saturation (or its damage). To meet this requirement, fast-gated (FG) single-photon avalanche diodes (SPADs) have been proposed in the literature [10,11]. The use of a FG-detector permits also the ability to reconstruct the DTOF with a larger dynamic range [12] and this great potentiality can be exploited if the laser source is powerful enough to increase the signal at late times, thus allowing the ability to probe deeper regions [13]. Additionally, the possibility to collect photons only at a given portion of the DTOF (by using a classical time-correlated single photon counting board or with simpler devices [14]) allows for the ability to increase the number of late photons, thus improving the visibility of possible deep inhomogeneities in the medium which affect the DTOF shape.

Different works based on FG-SPADs have been published, demonstrating experimentally that an improved spatial resolution, contrast and even chromophore quantification can be achieved using a gated detector [15–17]. However, to the best of our knowledge, the FG-SPAD detectors have never been adopted to demonstrate the improvement given by the use of short-SDD in fNIRS measurements. Indeed, the higher spatial resolution allows for the ability to retrieve more precisely the area wherein the activation occurs and the capability to collect photons in a given portion of the DTOF, improving the robustness of the detectability of a brain activation.

In this paper, we propose an instrument [18] based on a prototypal dual-channel (710 and 820 nm) high power (around 100 mW output for each channel) laser which allows for the ability to exploit at its best capacity the high-dynamic range acquisition technique. The laser light was spread over a sufficiently large area in order to guarantee skin safety by means of a spacer, while eye safety was guaranteed by means of the use of proper protection. Two state-of-the-art fast-gated SPAD modules (featuring a diffusion tail of about 90 ps [19]) were adopted to monitor both wavelengths, thus allowing for the recovery of the concentration of both O2Hb and HHB simultaneously. The proposed system was firstly characterized using some tests from two internationally-agreed protocols for performance assessment of TR-instrument/oximeter (such as the Instrument Response Function/IRF) from the "Basic Instrumental Protocol" (BIP) protocol [20], and contrast and contrast-to-noise ratio (CNR) from "nEUROPt" ones [21]). We then exploited the system for fNIRS measurements with the aim to demonstrate the improvement in terms of detectability of brain activation with respect to the classical oximetry instrument [22,23], simulating a free-running (i.e., not-gated) acquisition.

#### **2. Materials and Methods**

#### *2.1. Experimental Setup*

The setup used in the measurements is schematically depicted in Figure 1a. The laser source was a prototypal four-wave mixing laser (Fianium Ltd, Southampton, UK), providing two separate laser beams (710 and 820 nm) with a temporal duration of about 25 ps full width at half maximum (FWHM) at the repetition rate of 40 MHz. Each beam was properly attenuated by means of the variable optical attenuator (VOA, Edmund Optics Inc., Barrington, NJ, USA) and light was collimated into a 400 μm

fiber by means of a collimator. Collimated beams were then injected into the sample and retro-diffused photons were collected through two 1 mm core fibers (with a NA of 0.37) at a given source-detector distance (5 or 30 mm).

**Figure 1.** (**a**) Schematics of the setup used for the in vivo experiment. (**b**) Instrumental response function recorded for both fast-gated single-photon avalanche diode (SPAD) modules at their operating wavelengths (710 nm for fast-gated SPAD (FG-SPAD) 1—red line; and 820 nm for FG-SPAD 2—blue line). Time bin size: 3.05 ps. IRF = instrument response function, VOA = variable optical attenuator, TCSPC = time-correlated single photon counting. FWM: Four Wave Mixing.

In the in vivo measurements, it was fundamental to disentangle the contribution of the two wavelengths in order to recover the concentration of both oxygenated and deoxygenated hemoglobin. For this reason, an interferential filter (centered at 710 or 820 nm) was put at the other tip of each collection fiber. Such a solution allowed for the ability to separate the two beams and send each of them to a different fast-gated SPAD (FG-SPAD) module, thus having a parallel acquisition of the two wavelengths. The arrangement of the launching and detection fiber was conceived so as to probe the same area with both source-detector couples.

The FG-SPAD modules were state-of-the-art time-gated detectors, having an active area of 100 μm diameter and nearly 200 ps of raising edge for the gate opening. The full description of the FG-SPAD modules can be found in Boso et al. [10].

Concurrently with the optical pulse, the laser provided a synchronism signal, which was split into two paths. The former was fed to the FG-SPAD modules to trigger the enabling of the detector for a given period of time (5 ns). In order to acquire different "slices" of the DTOF, the detector was enabled at given delays with respect to the peak of the reflectance curve by means of a home-made delayer. It was a transmission-line based one, and it provided delays from 0 up to 6375 ps at a step of 25 ps.

On the other side, the second path of the split sync was sent to the time-correlated single photon counting (TCSPC) electronics, feeding the "stop" signal for both acquiring boards (SPC-130, Becker&Hickl GmbH, Berlin, Germany). The "start" signal for the TCSPC board was provided by each FG-SPAD module. Indeed, each module generated a pulse (timing OUT signal in Figure 1a) when a photon was detected within each SPAD, thus triggering an avalanche.

#### *2.2. Performance Assessment*

In order to impartially assess the performances of the proposed setup, we made use of some parameters of two internationally agreed protocols for diffuse optics instrumentation (the "Basic Instrumental Protocol" (BIP) [20] and "nEUROPt" [21] ones). Their aim is to characterize the TR setup for brain imaging.

For what concerns the BIP protocol, we acquired the IRF of the proposed setup at both wavelengths, coupling the launching and collection fiber with a diffusive medium (a thin film of TeflonTM) so as to excite all propagation modes of the detection fibers. In order to have a high-dynamic range measurement, we adopted the technique explained in Dalla Mora et al. and Tosi et al. [12,24], acquiring the re-emitted photons at different delays (from 0 up to 3450 ps at step of 25 ps) with an integration time of each delay of 10 s. Before rescaling the acquisitions and reconstructing the IRF, we subtracted the constant background (the so-called "dark count", which is not correlated with laser light) to all curves. We finally computed the FWHM and the dynamic range of the reconstructed IRF, as suggested by the BIP protocol.

We then characterized the proposed setup at 820 nm using two features of the "nEUROPt" protocol (contrast and contrast-to-noise ratio (CNR)) since they are meant to define quantitatively the capability of detecting a small absorption change within the probed volume (e.g., absorption changes in the cortex due to brain activation).

For both tests, we made use of a liquid phantom composed of a mix of water, Intralipid 20%, and Indian ink in such a proportion to have an absorption coefficient (μa) of 0.1 cm−1 and a reduced scattering one of 10 cm−1 at 820 nm. In order to simulate the absorption changes related to a localized brain activation, we made use of solid black cylindrical inclusions (height: 4.5 mm; diameter: 5.32 mm). In Martelli et al. [25], the equivalence relation between a realistic absorption change in a finite volume and the perturbation produced by a small black object had been demonstrated. In our case, we performed all tests using the 100 mm3 totally absorbing inclusion (which corresponded to a Δμa of 0.15 cm−1 over a 1 cm3 volume). The inclusion was posed beneath the surface of the phantom and its position was moved from 2 up to 30 mm in depth (the distance between the centroid of the inclusion and the surface). For each depth, the time-resolved curve was acquired in difference slices, delaying the opening of the gate from 0 ns (situation where the first re-emitted photon was recorded) up to 3250 ps at step of 250 ps. The VOA was set in order to have nearly 1 Mcps for each delay (except for the latest where it was not possible to increase the injected power). For each delay, we acquired five repetitions (1 s each) of the distribution of time-of-flight of re-emitted photons.

We computed contrast and CNR according the definition given in Wabnitz et al. [21] and briefly summarized in the following section.

The contrast is a primary assessment of the effect of variation in absorption coefficient. In the case of time-resolved reflectance measurements, often the analysis is performed using time-windows, meaning that the distribution of time-of-flight is subdivided into portions and contrast (as well as other features) and can be computed within those subdivisions of the curve. In case of the so-called "time-window" analysis, the contrast in each temporal window (*C*(*tW*)) was computed by the following equation:

$$\mathbf{C}^{(\mathbf{t\_w})} = \frac{\mathbf{N}\_0^{(\mathbf{t\_w})} - \mathbf{N}^{(\mathbf{t\_w})}}{\mathbf{N}\_0^{(\mathbf{t\_w})}},\tag{1}$$

where *N*(*tw*) 0 and *N*(*tw*) represent the overall number of counts computed in the time-window (*tw*) in the homogeneous or perturbed (e.g., when the inclusion mimicking the perturbation is inserted into the probed volume) state, respectively. In our case, each recorded curve was subdivided in time-windows of 500 ps in width. For the homogeneous acquisition, the inclusion was posed 30 mm aside from the source-detector couple, where its effect on the time-resolved curve was considered negligible.

The contrast-to-noise ratio affects the detection of a small variation in absorption coefficient. Indeed, CNR was defined in the case of "time-window" analysis as

$$\text{CNR}^{(\text{t}\_{\text{w}})} = \frac{\text{N}^{(\text{t}\_{\text{w}})} - \text{N}^{(\text{t}\_{\text{w}})}\_{0}}{\sigma \text{(N}^{(\text{t}\_{\text{w}})}\_{0})},\tag{2}$$

where σ *N*(*tw*) 0 is the fluctuation of counts in the unperturbed state (i.e., the inclusion far from the source-detector couple) computed over different repetitions. Usually, this standard deviation is mostly due to photonic noise, but also instrumental fluctuations (such as laser instability) can affect this value.

Both contrast and CNR were computed with and without the exploitation of the FG technique. Indeed, when analyzing the "non-gated" data, we used the measurements where the gate was opened before the first re-emitted photons were detected, and closed when the last ones were collected, thus mimicking the acquisitions obtained with a "classical" single-photon detector operating in the so-called "free-running" mode. On the other hand, for "gated" results, we analyzed the curve acquired enabling the FG-SPAD module with a given delay with respect to the peak of the reflectance curve.

Contrast and CNR were computed using the same portion of the curve (the same time-window defined with respect to the peak of the IRF). In this way the photons analyzed had reached the same mean depth (which was encoded in time for time-resolved measurements), thus leading to a fair comparison between the "gated" and "non-gated" technique.

#### *2.3. In Vivo Measurements*

Preliminary in vivo experiments were performed in order to test the performance of the system in a real case. Different finger-tapping experiments were performed to check system sensitivity in functional studies of the brain. Two healthy volunteers were recruited from the lab and informed consent was obtained (Subject one: male, 30 years old. Subject two: male, 43 years old). The competent ethics committee authorized the study.

The optical probe was placed over the left motor area and it was centered (according to the international 10–20 system for the Electro Encephalography -EEG- (electrode placement) at the C3 point. The measurement protocol was composed of 20 s of baseline, 20 s of motor task (finger opposition at a frequency rate as high as possible) and 20 s of recovery. This sequence was repeated five times with the right hand (controlateral hand) and five times with the left hand (ipsilateral hand). The entire protocol was repeated four times for each subject—using a short source-detector separation (ρ = 6 mm for both subjects) with and without the gating approach and using a long source-detector separation (ρ = 20 mm for subject 1 and ρ = 30 mm for subject 2) with and without the gating approach. The short source-detector separation distance was limited to 6 mm because of the size of the fiber terminator at the probe side, while the long source-detector separation was chosen in order to guarantee a photon

count as closer as possible to 1 Mcps. The same criterion was followed for the choice of the applied gating. Measurements were performed in a dimmed room to decrease the amount of background light.

In all the measurements, reemitted photons at both wavelengths were collected together with an integration time of 1 s. The photon time of flight distributions were divided into 25 time windows of 500 ps each. For each time window, we estimated the value of the quantity Δμavt at both the wavelengths [24], where Δμa is the variation of the absorption coefficient in the activated brain cortex, while v is the speed of light in the medium and t, the time spent by photons in the activated brain cortex. From this quantity at both wavelengths, we derived the changes in the amount of oxygenated hemoglobin (ΔO2Hb) and deoxygenated hemoglobin (ΔHHb), expressed in an arbitrary unit in the activated brain cortex. As mentioned in the phantom measurements section, we used the same time window for both gated and non-gated experiments in order to consider photons investigating the same depth. From the ΔO2Hb and ΔHHb curves, we can estimate contrast and CNR of the changes in the hemoglobin species with respect to the baseline.

#### **3. Results and Discussion**

#### *3.1. Characterization of the Setup (Phantom Measurements)*

Figure 1b represents the IRF of the proposed setup for the two channels. For both wavelengths the dynamic range was about 5.5 decades, while the computed FWHM was about 120 and 180 ps, respectively, for 710 and 820 nm. Some reflections were visible at 5.5 ns or 4.8 ns after the main peak (depending on the wavelength). Those reflections were most probably due to a second order of Fresnel reflections at the tip of the fibers. When conceiving the setup, the fibers' length was chosen in order to delay reflections at least some ns from the peak. Indeed, reflections did not represent a problem for the measurements, since they were quite late and several decades (about 4.5) under the main peak. Other small reflections presented at a shorter time were not significant since their amplitude was small enough to not distort the curve or saturate the dynamic range of the gated measurements.

Figure 2 reports the contrast and CNR computed at both the 5 and 30 mm source-detector distance (first and second row respectively) in the case of the fast-gating technique being applied ("gated", blue curve) or not ("non-gated", red lines). For all features, the time-window used for the analysis was opened 3 ns after the peak of the IRF (i.e., 3 ns after the injection of the light in the sample).

When a 5 mm source-detector distance was used, the fast-gating technique was fundamental to retrieve information about a buried perturbation. Indeed, the contrast recorded at 5 mm source-detector distance in the non-gated mode was low and very noisy, while, for the gated acquisitions, the contrast was much higher and its standard deviation (computed on the several repetitions) was nearly negligible. The improvement given by the fast-gating technique was clearly noticeable also in the CNR graph. Indeed, the curve relative to the gated acquisition was always higher than 10, even at the deepest position of the inclusion. On the other hand, the CNR for the non-gated acquisition was always less than 10 and, in most cases, also smaller than 1, meaning that the contrast was more sensitive to the noise of the measurements (e.g., photonic noise) than to the presence of an inclusion.

As predicted by the theory [9,26], the improvement given by the fast-gating technique was fundamental when using a small interfiber distance. Indeed, when the distance between the launching and detection fiber was shortened, a huge increase of the scarcely diffused photons (the so-called "early photons") was experienced. If a non-gated detector was used, this peak of "early photons" caused the saturation of its dynamic range, thus limiting the number of late photons detected and lowering the capability to discriminate a perturbation in depth.

The improvements given by the use of a fast-gated acquisition was less evident if a larger SDD was used, as it was clearly noticeable from curves reported in the last row of Figure 2. Indeed, the contrast obtained exploiting the fast-gating technique was slightly higher than what was achievable in the "non-gated" situation. The CNR graph clearly shows that in the gated mode there was a constant improvement of about a factor of 4.

**Figure 2.** Contrast (left column) and contrast-to-noise ratio (CNR) (right column) computed at two different source-detector distances (5 mm, top row; 30 mm, bottom row) with and without the exploitation of the fast-gating technique (blue and red curves, respectively).

This was mostly because, when using the fast-gating technique, we recorded photons enabling the FG-SPAD module with a delay of 2.25 ns with respect to what was done in the non-gated mode, thus rejecting the first part of the time-resolved curve. The delay was chosen in order to exploit the dynamic range of the detector with the latest possible delay. Since we collected the same number of photons (nearly 106 per acquisition) for both gated and not gated acquisition, it was obvious that in the former case there was a boost in the number of late photons detected, thus increasing the number of counts recorded in the region of the time-window used for analysis. This resulted in a decrease of the noise (which was also visible from the errorbar in the contrast graph) and in an increase of the CNR in the case of gated measurements.

Considering as a conservative threshold the 0.05 contrast and a CNR higher than 10, using the gated technique we could detect the inclusion until 2.6 cm and 3 cm in depth, with 5 and 30 mm source-detector distance, respectively.

This conclusion is in contrast with theoretical findings predicting a contrast that is always higher when using a small source-detector distance [26]. Such a disagreement can be due to both the effect of the "diffusion tail" of the detector and to the so-called "memory effect".

The former effect [19] is due to the non-perfectly efficient time-gating mechanism that causes a "leakage" of early photons in a temporal region of the curve where only the contribution of late photons are expected. This effect is as high as the source-detector distance decreases, since the peak of early photons increases. On the other side, the memory effect [27] causes a boost in the background level at longer delays that can hide the late photons, thus decreasing the maximum contrast achievable. Since the memory effect is proportional to the number of photons impinging on the detector when it is in the OFF state, its consequences are more evident when using a small source-detector distance.

#### *3.2. In Vivo Measurements*

Table 1 represents the summary of the in vivo experiments. Only data regarding the finger tapping with the contralateral hand were proposed because of the absence of any activation during finger tapping with the ipsilateral hand due to the lateralization of the motor area. For measurements at short SDD, contrast and CNR are higher for gated measurements, confirming that the time gated approach is necessary in order to perform time resolved functional measurements at short SDD. For long SDD measurements, this finding was not confirmed. In fact, only for subject 1 and for HHb in subject 2 did we have a slight improvement in contrast and CNR.

Figures 3 and 4 show the changes in hemoglobin species (Hb) concentrations averaged over the five repetitions of finger tapping with the contralateral hand for subject 1 and 2, respectively. Error bars represent the standard deviation over the five repetitions. Normally, brain activation is identified by a task related increase in O2Hb concentration and a corresponding decrease in HHb. We can clearly see this pattern in the gated column of both Figures 3 and 4, demonstrating the capability of a time gated system in performing functional measurements independently of the SDD, while in the non-gated approach this trend can be recognized only for long SDD.

**Figure 3.** In vivo measurements of subject 1. Changes in O2Hb (red) and HHb (blue) averaged over the five repetitions of finger tapping with the contralateral hand for the gated and non-gated approach (left and right column, respectively) at the short and long distance (upper and bottom lines, respectively). Error bars represent the standard deviation over the five repetitions.

**Table 1.** Summary of the in vivo measurements—finger tapping with the controlateral hand. Error for contrast and CNR is the standard deviation over the 5 finger tapping task repetitions. SDD = source-detector distance, O2Hb = oxygenated hemoglobin, HHb = deoxygenated hemoglobin.

**Figure 4.** In vivo measurements of subject 2. Changes in O2Hb (red) and HHb (blue) averaged over the five repetitions of finger tapping with the contralateral hand for the gated and non-gated approach (left and right column, respectively) at the short and long distance (upper and bottom lines, respectively). Error bars represent the standard deviation over the five repetitions.

#### **4. Conclusions**

The double-wavelength time-gated system proposed can perform measurements at SDD of few millimeters, showing the improvements theoretically predicted—increased contrast and CNR as highlighted by phantom and in vivo tests. Furthermore, the simultaneous acquisition of photon counts in a temporal window at two wavelengths allowed for the estimation of the changes in O2Hb and HHb concentration during brain activation. The use of time gating—necessary for short SDD—however, can be useful also for long SDD measurements. In fact, as demonstrated with phantom measurements, CNR for an absorption perturbation was always higher in time gated measurements. This was true because, in the presence of enough laser power, time gating allowed the exploitation of the entire dynamic range of the acquisition chain of the system to collect more late photons—which brought about information about the deeper structure in the medium—than without time gating. This last aspect was demonstrated with phantom measurements, but not completely with in vivo measurements. It is noteworthy that the in vivo measurements at different gates and at different distances were performed at different times; it is therefore plausible that their comparison is not completely reliable because of the intrinsic variability of the brain activation.

**Author Contributions:** Conceptualization, A.D.M., A.T., L.S., A.P. and D.C.; data curation, L.D.S., A.D.M. and D.C.; formal analysis, L.D.S., A.T., L.S., R.R. and D.C.; validation, L.D.S., A.D.M., A.T., L.S., A.P. and D.C.; writing—original draft, L.D.S., A.D.M., A.T., L.S., A.P. and D.C.; writing—review and editing, L.D.S., A.D.M., A.T., L.S., R.R., A.P. and D.C.

**Funding:** This research was partially supported from the EU project nEUROPt (grant agreement 201076).

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

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*


# **Geothermal Energy: Delivering on the Global Potential**

Edited by Paul L. Younger

Printed Edition of the Special Issue Published in *Energies*

www.mdpi.com/journal/energies

Paul L. Younger (Ed.)


**Abstract**

The aim of this study is to analyze the limitations of the Bilingual Program of the Community of Madrid and present proposals for improvement to solve them. The methodology used is based on the collection of information through different sources. The limitations indicated are as follows: the attention to the students with specific need of educational support; the coordination of teachers; teacher training; English language assistants; human resources; material resources; the motivation of students, teachers, and English language assistants; the methodology; the loss of content and vocabulary in the mother tongue; the need for language academies or extracurricular activities in English to complete training; and segregation. It is concluded that the greatest and most urgent limitation of the program is the attention of students with special educational needs. Proposals for improvement are presented, especially in relation to training and coordination with specialist teachers.

**Keywords:** bilingual education, bilingual program, limitations, proposals for improvement, primary education


**3. Measurement system and dedicated hardware**

A measurement system has been built in order to verify the approach. Particular components, like the UWB-antennas used and the working principle of the whole set-up, are discussed in this section. A dedicated time domain transmission oscilloscope is presented. It is capable of transmitting and receiving UWB pulses with several gigahertz of bandwidth. It is specially tailored to the measurement problem and therefore functions with less hardware and software, is far more compact, and is cheaper compared to a universal laboratory instrument.

### **3.1. Compact ultra-wideband antennas and arrays**

The sensor system should be able to transmit and receive ultra-wideband (UWB) signals in two orthogonal polarizations. This is of great importance if the orientation of the object

4 Will-be-set-by-IN-TECH 326 Ultra-Wideband Radio Technologies for Communications, Localization and Sensor Applications ISOPerm: Non-Contacting Measurement of Dielectric Properties of Irregular Shaped Objects 5

(a) Dual polarised dipole. (b) Dual polarised array.

antenna matching. The radiating elements have a maximum diameter of 30mm. The feed is provided through a coaxial-slotline transition, so that no other components are necessary. The distance between two dipoles having an equal polarization is 80mm in each direction; between dipoles with the same polarization it is √2 × 80mm. The dimensions are overall 160mm × 320mm × 100mm. The matching is better than 10dB between 1.6GHz and 4.2GHz.

ISOPerm: Non-Contacting Measurement of Dielectric Properties of Irregular Shaped Objects 327

The core of the measurement system is the signal generation and sampling. The requirements for the system are that the step signals generated have a rise time in the range of 100ps, and a large amplitude. These signals are transmitted through the measurement path after which, the step response has to be sampled. The system should have a high bandwidth, low noise and low jitter. Furthermore it should be as compact and affordable as possible. Therefore, an impulse technique using equivalent time sampling is the method of choice. Classical swept sine wave techniques, as used in network analysers or real time digital sampling oscilloscopes,

Equivalent time sampling is based on the repetitive stimulation by a measurement signal with a cycle duration of *T*0. At every repetition of the measurement signal the moment of sampling is shifted by Δ*T*. The cycle of the sampling clock is then *T*1 = *T*0 + Δ*T*. The

conversion and the data transport and storage is greatly reduced compared to a real time oscilloscope. This technique has been used for decades; with recently available MMICs, both high performance and cheap hardware can be achieved. It is the method of choice for the sampling of signals with huge bandwidth in the GHz range with a high resolution. Laboratory instruments employing this technique have been used for preliminary investigations. In order to demonstrate the system performance and accuracy under practical conditions, dedicated hardware was developed. Two specially tailored systems were investigated. They differed

Δ*T T*1. The effort regarding analog-to-digital

The crosstalk between the single elements is maximum −20dB at 10GHz.

**3.2. Signal generation and sampling**

are too expensive, complex, and bulky.1

track and hold amplifier

XTAL

HP

f2 f2

output generator

reference path

input

sampled signal is therefore stretched to *T*2 = *T*0

ADC FIFO

f2

clock generation

step

f1

memory

microcontroller

computer

**Figure 3.** Time domain transmission measurement system with a microcontroller.

measurements carried out with this method which lead to excellent results [20].

(a) Block diagram. (b) Prototype: integrated on one PCB except microcontroller and computer.

1 A rather exotic concept employing M-sequences can be found in this book (see chapter HALOS). There are also

under test (e.g. prolate ellipsoid) is not known. There are several polarization set-ups for receiving co- and cross-polarized parts of the scattered wave. In the set-up used latterly the polarization of the transmitting antenna is tilted at 45◦ with respect to the receiving array so that ideally, without an object under test, the receiving antennas receive the transmitted pulse in equal proportions for each polarization. Therefore, dual orthogonal UWB antennas have been developed. Furthermore the single antennas should be as compact as possible while operating in the lower gigahertz range. The arrays presented later are not group antennas as commonly used. The received signal of every antenna is sampled separately in order not to confound the information at different locations with respect to the object under test. Operating at frequencies in the lower GHz range is anticipated because the electromagnetic waves have a higher penetration depth (especially for lossy dielectrics with a high permittivity) and therefore will be more affected.

There are many UWB antennas suitable for building dual polarized antennas [21–25]. Dual polarized dipole (see [26–31]) and horn antennas (see [32–35]) can already be found in the literature. Most of these antennas have the disadvantage that they have crossed feed points. It is possible to overcome this problem with the configuration depicted in Figure 2. Both the four radiation elements and the whole configuration are two-fold symmetric. Two dipoles are excited in even mode in horizontal as well as in vertical polarization. Ideally, due to the symmetry, both planes are decoupled. The diameter of the radiation elements is 24mm; the whole PCB containing the four radiating elements measures 50mm × 50mm. The antenna is equipped with two orthogonal feed networks; each is composed of a two stage Wilkinson divider and two tapered baluns. The horizontal and vertical polarization feed networks have dimensions of 50mm × 142mm and 50mm × 103mm, respectively. In order to avoid backward radiation, absorber material is mounted appropriately (the absorber material was removed for the photograph shown in Figure 2).

The matching at both ports is better than 7dB in a frequency range from 2GHz to 5.7GHz; the isolation is better than 30dB. At 5GHz the radiated crosspolar portion is suppressed by more than 25dB (boresight); the gain is 5.2dB. Measurements with a line containing four of these antennas lead to excellent results but the whole antenna structure is rather complex. Since a scenario with an object moving orthogonally to the line array is anticipated (e.g. a conveyor belt), a simpler arrangement can be used. The antenna array shown in Figure 2 consists of eight single linear broadband dipoles. The geometry is optimized regarding the antenna matching. The radiating elements have a maximum diameter of 30mm. The feed is provided through a coaxial-slotline transition, so that no other components are necessary. The distance between two dipoles having an equal polarization is 80mm in each direction; between dipoles with the same polarization it is √2 × 80mm. The dimensions are overall 160mm × 320mm × 100mm. The matching is better than 10dB between 1.6GHz and 4.2GHz. The crosstalk between the single elements is maximum −20dB at 10GHz.

### **3.2. Signal generation and sampling**

4 Will-be-set-by-IN-TECH

(a) Dual polarised dipole. (b) Dual polarised array.

under test (e.g. prolate ellipsoid) is not known. There are several polarization set-ups for receiving co- and cross-polarized parts of the scattered wave. In the set-up used latterly the polarization of the transmitting antenna is tilted at 45◦ with respect to the receiving array so that ideally, without an object under test, the receiving antennas receive the transmitted pulse in equal proportions for each polarization. Therefore, dual orthogonal UWB antennas have been developed. Furthermore the single antennas should be as compact as possible while operating in the lower gigahertz range. The arrays presented later are not group antennas as commonly used. The received signal of every antenna is sampled separately in order not to confound the information at different locations with respect to the object under test. Operating at frequencies in the lower GHz range is anticipated because the electromagnetic waves have a higher penetration depth (especially for lossy dielectrics with a high permittivity) and

There are many UWB antennas suitable for building dual polarized antennas [21–25]. Dual polarized dipole (see [26–31]) and horn antennas (see [32–35]) can already be found in the literature. Most of these antennas have the disadvantage that they have crossed feed points. It is possible to overcome this problem with the configuration depicted in Figure 2. Both the four radiation elements and the whole configuration are two-fold symmetric. Two dipoles are excited in even mode in horizontal as well as in vertical polarization. Ideally, due to the symmetry, both planes are decoupled. The diameter of the radiation elements is 24mm; the whole PCB containing the four radiating elements measures 50mm × 50mm. The antenna is equipped with two orthogonal feed networks; each is composed of a two stage Wilkinson divider and two tapered baluns. The horizontal and vertical polarization feed networks have dimensions of 50mm × 142mm and 50mm × 103mm, respectively. In order to avoid backward radiation, absorber material is mounted appropriately (the absorber material was removed for

The matching at both ports is better than 7dB in a frequency range from 2GHz to 5.7GHz; the isolation is better than 30dB. At 5GHz the radiated crosspolar portion is suppressed by more than 25dB (boresight); the gain is 5.2dB. Measurements with a line containing four of these antennas lead to excellent results but the whole antenna structure is rather complex. Since a scenario with an object moving orthogonally to the line array is anticipated (e.g. a conveyor belt), a simpler arrangement can be used. The antenna array shown in Figure 2 consists of eight single linear broadband dipoles. The geometry is optimized regarding the

**Figure 2.** Ultra-wideband antenna structures.

therefore will be more affected.

the photograph shown in Figure 2).

The core of the measurement system is the signal generation and sampling. The requirements for the system are that the step signals generated have a rise time in the range of 100ps, and a large amplitude. These signals are transmitted through the measurement path after which, the step response has to be sampled. The system should have a high bandwidth, low noise and low jitter. Furthermore it should be as compact and affordable as possible. Therefore, an impulse technique using equivalent time sampling is the method of choice. Classical swept sine wave techniques, as used in network analysers or real time digital sampling oscilloscopes, are too expensive, complex, and bulky.1

Equivalent time sampling is based on the repetitive stimulation by a measurement signal with a cycle duration of *T*0. At every repetition of the measurement signal the moment of sampling is shifted by Δ*T*. The cycle of the sampling clock is then *T*1 = *T*0 + Δ*T*. The sampled signal is therefore stretched to *T*2 = *T*0 Δ*T T*1. The effort regarding analog-to-digital conversion and the data transport and storage is greatly reduced compared to a real time oscilloscope. This technique has been used for decades; with recently available MMICs, both high performance and cheap hardware can be achieved. It is the method of choice for the sampling of signals with huge bandwidth in the GHz range with a high resolution. Laboratory instruments employing this technique have been used for preliminary investigations. In order to demonstrate the system performance and accuracy under practical conditions, dedicated hardware was developed. Two specially tailored systems were investigated. They differed

and computer.

**Figure 3.** Time domain transmission measurement system with a microcontroller.

1 A rather exotic concept employing M-sequences can be found in this book (see chapter HALOS). There are also measurements carried out with this method which lead to excellent results [20].

0 100 200 300 400 500

TDS8000 TDT

0 100 200 300 400 500

TDS8000 TDT

t [ps]

(b) Monocycle.

−1

−0.5

0

normalized amplitude

**Figure 5.** Comparison between the proposed dedicated hardware (TDT) and a Textronix TDS8000 (no

TDT

(a) Block diagram. (b) Photograph.

of Δ*T* = 4ps. For the measurements presented later a much lower resolution is sufficient, e.g.

The entire measurement system is shown is Figure 6. The step signal of the TDT is amplified by a power amplifier (PA) before being transmitted through a horn antenna. The object under test is illuminated by the transmitted signal and portions of the scattered signal are received by an array of eight receiving antennas, which are arranged in two orthogonal polarizations. In order to emulate a 2-dimensional array, the movable table is moved orthogonally to the orientation of the receiving array. It is crucial to the multivariate calibration, which is applied

LNA

PA

combiner or switch

0.5

1

ISOPerm: Non-Contacting Measurement of Dielectric Properties of Irregular Shaped Objects 329

t [ps]

(a) Pulse.

movable table

0 0.2 0.4 0.6 0.8 1

normalized amplitude

averaging).

object under test

Tx

Rx

Δ*T* = 40ps.

**Figure 6.** Entire measurement set-up.

**3.3. Entire measurement setup**

**Figure 4.** Time domain transmission measurement system with a FPGA.

mostly in the digital part behind the analog-digital converter (ADC); the front end was very similar.

A block diagram and a photograph of the first system using a microcontroller are shown in Figure 3. Two slightly detuned clocks are synthesized from a crystal. One clock triggers a step generator, which transmits a step signal with about 30ps rise time (20% to 80%) and up to 3V amplitude. The other triggers the track and hold amplifier (having an input bandwidth of 13GHz), the 12bit ADC, and the asynchronous FIFO (first in first out) memory with a width of 18bit and a depth of 32k. A microcontroller and a computer control the system. The FIFO memory can be read out, and the clock generation can be programmed via the I2C-bus. Clock signals up to hundreds of megahertz can be chosen with an accuracy of 1ppm. The dimensions of the RF-PCB are 90mm × 60mm. A reference path is necessary to measure time differences, because the phase of the two clock signals is not captured. A FIFO has to be used because the microcontroller is not able to read in 12bit words at tens of megahertz. Therefore, one waveform of up to 32k is captured at a time and is then transferred to the computer. Averaging, in order to improve the SNR, is carried out on the computer.

An improvement is the application of a field programmable gate array (FPGA). A modular system is shown in Figure 4. There is no more need for a fast external FIFO because data can be read into the FPGA directly up to about 66MHz. Furthermore, no reference path is needed, because the FPGA is able to discriminate the phase between both clock signals. Up to 212 = 4096 sampling points are possible while an averaging of up to 213 = 8192 waveforms can be carried out on the FPGA (this is restricted due to the internal resources of the FPGA used2).

Since both front ends are similar they offer a comparable performance. The RMS noise level of the receiver is < 0.9mV, the RMS jitter is < 0.7ps (no averaging). Figure 5 shows a comparison of sampled test signals between the proposed hardware and a Tektronix TDS8000. The frequencies are chosen as *f*1 = 50.01MHz and *f*2 = 50MHz, which leads to a resolution

2 A Spartan 6 from Xilinx is used.

**Figure 5.** Comparison between the proposed dedicated hardware (TDT) and a Textronix TDS8000 (no averaging).

6 Will-be-set-by-IN-TECH

FPGA

computer

mostly in the digital part behind the analog-digital converter (ADC); the front end was very

A block diagram and a photograph of the first system using a microcontroller are shown in Figure 3. Two slightly detuned clocks are synthesized from a crystal. One clock triggers a step generator, which transmits a step signal with about 30ps rise time (20% to 80%) and up to 3V amplitude. The other triggers the track and hold amplifier (having an input bandwidth of 13GHz), the 12bit ADC, and the asynchronous FIFO (first in first out) memory with a width of 18bit and a depth of 32k. A microcontroller and a computer control the system. The FIFO memory can be read out, and the clock generation can be programmed via the I2C-bus. Clock signals up to hundreds of megahertz can be chosen with an accuracy of 1ppm. The dimensions of the RF-PCB are 90mm × 60mm. A reference path is necessary to measure time differences, because the phase of the two clock signals is not captured. A FIFO has to be used because the microcontroller is not able to read in 12bit words at tens of megahertz. Therefore, one waveform of up to 32k is captured at a time and is then transferred to the computer.

An improvement is the application of a field programmable gate array (FPGA). A modular system is shown in Figure 4. There is no more need for a fast external FIFO because data can be read into the FPGA directly up to about 66MHz. Furthermore, no reference path is needed, because the FPGA is able to discriminate the phase between both clock signals. Up to 212 = 4096 sampling points are possible while an averaging of up to 213 = 8192 waveforms can be carried out on the FPGA (this is restricted due to the internal resources of the FPGA

Since both front ends are similar they offer a comparable performance. The RMS noise level of the receiver is < 0.9mV, the RMS jitter is < 0.7ps (no averaging). Figure 5 shows a comparison of sampled test signals between the proposed hardware and a Tektronix TDS8000. The frequencies are chosen as *f*1 = 50.01MHz and *f*2 = 50MHz, which leads to a resolution

(a) Block diagram. (b) Prototype: modular system except computer.

track and hold amplifier

XTAL

f2 f2

output generator

ADC

clock generation f1 f2

step

**Figure 4.** Time domain transmission measurement system with a FPGA.

Averaging, in order to improve the SNR, is carried out on the computer.

f1

input

similar.

used2).

2 A Spartan 6 from Xilinx is used.

of Δ*T* = 4ps. For the measurements presented later a much lower resolution is sufficient, e.g. Δ*T* = 40ps.

### **3.3. Entire measurement setup**

The entire measurement system is shown is Figure 6. The step signal of the TDT is amplified by a power amplifier (PA) before being transmitted through a horn antenna. The object under test is illuminated by the transmitted signal and portions of the scattered signal are received by an array of eight receiving antennas, which are arranged in two orthogonal polarizations. In order to emulate a 2-dimensional array, the movable table is moved orthogonally to the orientation of the receiving array. It is crucial to the multivariate calibration, which is applied

#### 8 Will-be-set-by-IN-TECH 330 Ultra-Wideband Radio Technologies for Communications, Localization and Sensor Applications ISOPerm: Non-Contacting Measurement of Dielectric Properties of Irregular Shaped Objects 9

later, that the received signals of each antenna placed at its unique position are sampled separately. Therefore, the antennas are switched sequentially onto the input of the TDT and are sampled there, after being amplified by a low noise amplifier (LNA). The MOSFET switch used requires relatively high hard- and software effort and it has to be synchronised with the other components included in the system. Furthermore, these switches have a latency of some nanoseconds. The switch has an insertion loss of 4dB at 2GHz and 8.6dB at 8GHz.

A simpler solution is to use a broadband combiner. The signals received by the individual antennas are combined using a broadband eight-way Wilkinson divider. Prior to combining there is a delay of *τ* = 2ns between two adjacent inputs. It is possible to separate the individual pulses in time. Compared to using a receiver having more channels, or using a switch, the hardware effort is greatly reduced and the instantaneous sampling of all pulses is possible in about 20ns. One of the disadvantages that has to be considered is that there is an increased insertion loss of ideally 9dB; at 1GHz and 5GHz the measured insertion losses are 9.3dB and 10.7dB, respectively. There are also losses and temperature dependencies due to the delay lines, which have to be taken into account.

(a) One of the moulds manufactured from

40 50 60

4.88 % moisture 20.93 % moisture

**Figure 8.** Examples of received pulses from two objects having different moisture content.

multivariate calibration methods which will be explained in the next section 5.

of liquid, and another ten with about 140g. Therefore, overall 20 objects were measured. A photograph of one of the measured objects containing 113.15g ethanol and 27.94g of water is

The time domain signals are gated in order to extract the time interval of interest: about 10 to 30 equidistant amplitude values per pulse are empirically chosen and subjected to the

The sequence of the eight sampled pulses consists of 217 sampling points. There are four positions of the movable table which then leads to 4 × 217 = 868 data points per object. These pulses are modified by many factors: the shape, the position, the rotation, and the intrinsic variables of the material under test. However, as can be seen in Figure 8 the value to be measured has only a relatively low influence on the apparent shape of the pulses, and

−0.15 −0.1 −0.05 0 0.05 0.1 0.15

amplitude [V]

t [ns]

(a) Instantaneous sampled pulses received by the eight antennas using a broadband Wilkinson divider and

**5. Multivariate data processing**

(b) Mould filled with moist clay granules. (c) A bottle filled

ISOPerm: Non-Contacting Measurement of Dielectric Properties of Irregular Shaped Objects 331

with an ethanol-water

mixture.

40.5 41 41.5

t [ns]

(b) Closer look at the first received pulse.

Styrofoam.

−0.2

delay lines.

shown in Figure 7.

−0.1

amplitude [V]

0

0.1

0.2

**Figure 7.** Objects under test.


**3. Multi-Agent Simulator**

This section briefly describes the basics of multi-agent simulation (MAS). NetLOGO was chosen as the tool for the implementation of the MAS, which will be explained further below. In Section 3.3, the most important elements of the programmed MAS are presented concisely.


**6.2 Familial linkage**

In the Collaborative Longitudinal Evaluation of Keratoconus (CLEK) study in a cohort of 1209 patients mixed in ethnicity in the USA, 829 (69%) white, 240 (20%) black, and 99 (8%) Hispanic, family history of keratoconus was reported in 13.5% of the cases [65]. After follow-up for 8 years, the inheritance patterns were not established [66]. In another study from the USA, more females reported family history than males. But it was unclear whether it was a difference in attitude on reporting or a genuine gender difference in familial link [57]. In the New Zealand study, familial aggregation analysis showed keratoconus familial rate of 23.5% [59]. In Scotland, family history occurred in 5% of 186 white patients [67]. In a report from North India among 120 patients, 6 (5%) had family history [68]. In a review of keratoconus in Asians, family history ranged from 4.4 to 23.5% [40]. Overall, reported family history of keratoconus has widely ranged from 6% up to about 25% [37, 65].

A recent systemic review and meta-analysis on 29 eligible reports from different parts of the world selected from 3996 articles revealed family history as the strongest risk factor (odds ratio 6.42; 95% CI:2.59–10.24) among other established risk factors: eye rubbing (odds ratio 3.09; 95% CI:2.17–4.00), eczema (odds ratio 2.95; 95% CI:1.30–4.59), asthma (odds ratio 1.94; 95% CI:1.30–2.58), and allergy (odds ratio 1.42; 95% CI:1.06–1.79) [69]. Overall, no obvious differences exist in clinical or ophthalmic presentations between familial and sporadic keratoconus. As a genetic disease, keratoconus has shown weak penetrance with variable expressions. While family aggregation and linkage studies showed familial inheritance, the majority of reported keratoconus patients are sporadic.

#### **7. Mapping the keratoconus genes**

Candidate gene analysis, family linkage analysis, and more recently genomewide association study with support by candidate gene association studies and next generation sequencing, have contributed to identification of genetic loci or gene variants that are in association with keratoconus [70]. All are genetic associations. No direct keratoconus causing mutation has been identified.

#### **7.1 Linkage analysis**

Single nucleotide polymorphisms (SNPs) and microsatellite markers covering the whole genome have been attempted to find keratoconus loci or even genes in family pedigrees and sib pairs. But even in families, penetrance of keratoconus is low and clinical presentations are variable. A large number of samples have to be available. Vigorous research among various ethnic groups have identified a number of keratoconus loci which are replicable and of maximum LOD score greater than 3, in a collection of 67 sib pair Hispanic families, two-stage genome-wide analysis of 380 microsatellite markers in totally 351 subjects, among them 110 were affected by keratoconus which has led to identification of a multitude of loci in chromosomes *2q, 3p, 4q, 5q31, 5p, 9p, 9q34, 11p, 12p, 14q*, and *17q* [71]. In a collection of 25 Italian families, genome-wide scan of microsatellite markers in 77 affected and 57 unaffected members have mapped chromosomal regions for keratoconus in *5q32-q33; 5q21.2, 14q11.2,* and *15q2.32* [72]. Some of these loci had been replicated or refined in further investigations: *2p24* [73]; *3p14-q13* [74]; *5q14.3-q21.1* [75]; 5q31.1-q35.3 [72]; *13q32* [76]; *16q22.3-q23.1* [77]; and *20q12* [78]. These studies were mostly on European populations. The large number of associated loci exemplified the genetic heterogeneity of keratoconus. No disease-causing mutation has been identified from these loci.

#### **7.2 Genome-wide association studies (GWAS)**

Most GWAS in connection with keratoconus were conducted on the two quantitative traits of central corneal thickness (CCT) and corneal curvature (CC). These two are characteristic, but not exclusive, traits of keratoconus, with the corresponding morphological features of central corneal thinning and corneal steepness.

#### *7.2.1 Central corneal thickness (CCT)*

There were nine reported GWAS on central corneal thickness (CCT) up to July 2017. They were summarized in an excellent review [70]. All involved meta-analysis

**61**

*7.2.3 Other approaches*

*Molecular Genetics of Keratoconus: Clinical Implications DOI: http://dx.doi.org/10.5772/intechopen.90623*

the highest among all reported SNPs [79].

and 4.0 × 10−5

P-values of 2.5 × 10−7

*7.2.2 Corneal curvature (CC)*

ethnicities [86].

in initial and validation cohorts. More than 40 SNPs in 30 genes were reported. The biggest sample size was with totally 13,057 European and 6963 Asian samples, while the primary cohort was comprised of 874 patients and 6085 controls. A multitude of keratoconus associated SNPs in 26 loci was identified. SNP rs9938149 in the *BANP-ZNF469* gene attained GWAS significance of P-value = 2.4 × 10−49, which was

Some of the gene variants were analyzed in keratoconus and control cohorts. Possibly there was ethnic specificity. Rs96067 in COL8A2 was identified in a Singaporean cohort of 2538 Indians and 2542 Malays [80] but not in a separate GWAS study 0f European study subjects with 3931 German and 1418 Dutch study subjects. However, SNPs in *BANP-ZNF469* and COL5A1 were replicated [81]. In a cohort of Australian white study subjects with 933 keratoconus patients and more than 4000 controls, keratoconus susceptibility was detected at the *HGF* locus [82]. The risk factor allele rs3735520 was associated with keratoconus in a Czech cohort of 165 patients and 193 controls [83] and Australian whites of 157 patients and 673 controls [84]. Another study involved two independent cohorts of keratoconus patients, involving 222 Caucasian patients, 687 African Americans, 3324 Caucasian controls and 307 individuals from 70 keratoconus families reported strong association of *Lysyl Oxidase* gene (*LOX)* polymorphisms with keratoconus, with meta

respectively [85]. In a meta-analysis of 14 studies comprising 17,803 individuals of European ancestry 44 loci associated with CCT were found, two of them were LTBP1 and WNT10A 42 European specific while the rest occurred also in Asian

that contributed to identification of susceptibility genes. Eight SNPs in MTOR/ FRAP1 and 7 SNPs in PDGFRA were found in 10,008 samples from three population groups in Singapore: Malays, Chinese and Indians [87]. Another big cohort of 12,660 Asians included Japanese in addition to Malays, Chinese, and Indians, [88] Associations of CMPK1 rs17103186 and RBP3 rs11204213 with CC, and RBP3 rs11204213 with axial length were discovered. In the Avon Longitudinal Study of Parents and Children (ALSPAC) cohort of 2023 individuals of white European descent, rs6554163 in *PDGFRA* was associated with both CC and axial length [89]. In a cohort of 1013 Australian individuals, 1788 twins and their families, CC was associate with rs2114039 near *PDGFRA* and rs2444240 which is 31 kb upstream to *TRIM29*, but not with any SNP at the FRAP1 locus [90]. In totally 15,168 study samples that included Japanese, Chinese or European ethicity, rs10453441 in WNT7B was strongly associated with both CC and axial length [91]. WNT7B rs10453441 is the only SNP associated with both traits of CCT and CC [92]. While these SNPs have no consistent and strong association with keratoconus, an exome sequencing analysis identified a *WNT10A* variant that was associated with corneal thickness and keratoconus [93]. In contrast, a GWAS on direct association of

Six GWAS on corneal curvature (CC) have been reported in multi-ethnic cohorts

keratoconus patients and controls involving 222 patients and 3324 controls found no

Apart from corneal morphologic features, a recent GWAS investigated corneal biomechanical properties in an European cohort of 6645 participants and 2384 participants from the British TwinsUK study [36]. The second stage of the association

GWAS significance of associated gene variants [94].

for LOX SNPs rs2956540 and rs10519694

*Ocular Surface Diseases - Some Current Date on Tear Film Problem and Keratoconic Diagnosis*

in Asians, family history ranged from 4.4 to 23.5% [40]. Overall, reported family history of keratoconus has widely ranged from 6% up to about 25% [37, 65].

reported keratoconus patients are sporadic.

No direct keratoconus causing mutation has been identified.

**7.2 Genome-wide association studies (GWAS)**

*7.2.1 Central corneal thickness (CCT)*

**7. Mapping the keratoconus genes**

**7.1 Linkage analysis**

A recent systemic review and meta-analysis on 29 eligible reports from different parts of the world selected from 3996 articles revealed family history as the strongest risk factor (odds ratio 6.42; 95% CI:2.59–10.24) among other established risk factors: eye rubbing (odds ratio 3.09; 95% CI:2.17–4.00), eczema (odds ratio 2.95; 95% CI:1.30–4.59), asthma (odds ratio 1.94; 95% CI:1.30–2.58), and allergy (odds ratio 1.42; 95% CI:1.06–1.79) [69]. Overall, no obvious differences exist in clinical or ophthalmic presentations between familial and sporadic keratoconus. As a genetic disease, keratoconus has shown weak penetrance with variable expressions. While family aggregation and linkage studies showed familial inheritance, the majority of

Candidate gene analysis, family linkage analysis, and more recently genomewide association study with support by candidate gene association studies and next generation sequencing, have contributed to identification of genetic loci or gene variants that are in association with keratoconus [70]. All are genetic associations.

Single nucleotide polymorphisms (SNPs) and microsatellite markers covering the whole genome have been attempted to find keratoconus loci or even genes in family pedigrees and sib pairs. But even in families, penetrance of keratoconus is low and clinical presentations are variable. A large number of samples have to be available. Vigorous research among various ethnic groups have identified a number of keratoconus loci which are replicable and of maximum LOD score greater than 3, in a collection of 67 sib pair Hispanic families, two-stage genome-wide analysis of 380 microsatellite markers in totally 351 subjects, among them 110 were affected by keratoconus which has led to identification of a multitude of loci in chromosomes *2q, 3p, 4q, 5q31, 5p, 9p, 9q34, 11p, 12p, 14q*, and *17q* [71]. In a collection of 25 Italian families, genome-wide scan of microsatellite markers in 77 affected and 57 unaffected members have mapped chromosomal regions for keratoconus in *5q32-q33; 5q21.2, 14q11.2,* and *15q2.32* [72]. Some of these loci had been replicated or refined in further investigations: *2p24* [73]; *3p14-q13* [74]; *5q14.3-q21.1* [75]; 5q31.1-q35.3 [72]; *13q32* [76]; *16q22.3-q23.1* [77]; and *20q12* [78]. These studies were mostly on European populations. The large number of associated loci exemplified the genetic heterogeneity of keratoconus. No disease-causing mutation has been identified from these loci.

Most GWAS in connection with keratoconus were conducted on the two quantitative traits of central corneal thickness (CCT) and corneal curvature (CC). These two are characteristic, but not exclusive, traits of keratoconus, with the corresponding morphological features of central corneal thinning and corneal steepness.

There were nine reported GWAS on central corneal thickness (CCT) up to July 2017. They were summarized in an excellent review [70]. All involved meta-analysis

**60**

in initial and validation cohorts. More than 40 SNPs in 30 genes were reported. The biggest sample size was with totally 13,057 European and 6963 Asian samples, while the primary cohort was comprised of 874 patients and 6085 controls. A multitude of keratoconus associated SNPs in 26 loci was identified. SNP rs9938149 in the *BANP-ZNF469* gene attained GWAS significance of P-value = 2.4 × 10−49, which was the highest among all reported SNPs [79].

Some of the gene variants were analyzed in keratoconus and control cohorts. Possibly there was ethnic specificity. Rs96067 in COL8A2 was identified in a Singaporean cohort of 2538 Indians and 2542 Malays [80] but not in a separate GWAS study 0f European study subjects with 3931 German and 1418 Dutch study subjects. However, SNPs in *BANP-ZNF469* and COL5A1 were replicated [81]. In a cohort of Australian white study subjects with 933 keratoconus patients and more than 4000 controls, keratoconus susceptibility was detected at the *HGF* locus [82]. The risk factor allele rs3735520 was associated with keratoconus in a Czech cohort of 165 patients and 193 controls [83] and Australian whites of 157 patients and 673 controls [84]. Another study involved two independent cohorts of keratoconus patients, involving 222 Caucasian patients, 687 African Americans, 3324 Caucasian controls and 307 individuals from 70 keratoconus families reported strong association of *Lysyl Oxidase* gene (*LOX)* polymorphisms with keratoconus, with meta P-values of 2.5 × 10−7 and 4.0 × 10−5 for LOX SNPs rs2956540 and rs10519694 respectively [85]. In a meta-analysis of 14 studies comprising 17,803 individuals of European ancestry 44 loci associated with CCT were found, two of them were LTBP1 and WNT10A 42 European specific while the rest occurred also in Asian ethnicities [86].

### *7.2.2 Corneal curvature (CC)*

Six GWAS on corneal curvature (CC) have been reported in multi-ethnic cohorts that contributed to identification of susceptibility genes. Eight SNPs in MTOR/ FRAP1 and 7 SNPs in PDGFRA were found in 10,008 samples from three population groups in Singapore: Malays, Chinese and Indians [87]. Another big cohort of 12,660 Asians included Japanese in addition to Malays, Chinese, and Indians, [88] Associations of CMPK1 rs17103186 and RBP3 rs11204213 with CC, and RBP3 rs11204213 with axial length were discovered. In the Avon Longitudinal Study of Parents and Children (ALSPAC) cohort of 2023 individuals of white European descent, rs6554163 in *PDGFRA* was associated with both CC and axial length [89]. In a cohort of 1013 Australian individuals, 1788 twins and their families, CC was associate with rs2114039 near *PDGFRA* and rs2444240 which is 31 kb upstream to *TRIM29*, but not with any SNP at the FRAP1 locus [90]. In totally 15,168 study samples that included Japanese, Chinese or European ethicity, rs10453441 in WNT7B was strongly associated with both CC and axial length [91]. WNT7B rs10453441 is the only SNP associated with both traits of CCT and CC [92]. While these SNPs have no consistent and strong association with keratoconus, an exome sequencing analysis identified a *WNT10A* variant that was associated with corneal thickness and keratoconus [93]. In contrast, a GWAS on direct association of keratoconus patients and controls involving 222 patients and 3324 controls found no GWAS significance of associated gene variants [94].

### *7.2.3 Other approaches*

Apart from corneal morphologic features, a recent GWAS investigated corneal biomechanical properties in an European cohort of 6645 participants and 2384 participants from the British TwinsUK study [36]. The second stage of the association

study involved 752 keratoconus patients as compared with 974 British TwinsUK or 13,828 EPIC-Norfolk. The results showed a likely role in development of keratoconus with 5 associated loci in CH, *ANAPC1*, *ADAMTS8*, *ADAMTS17*, *ABCA6*, and *COL6A1* [36].

It is notable that there are a lot more keratoconus genes that are identified through studies on quantitative traits, especially central cornea thickness and cornea curvature, than on direct association with keratoconus against controls (**Figure 1**).


**1. Introduction**

In recent years, heavy-duty vehicle transportation has become important all around the world [1,2]. The concept of dynamic interactions between the vehicle and the road surface has been widely studied, and a series of problems caused by the coupling vibrations of the vehicle structure and road have become very prominent. Thanks to the in-depth study of this issue by scholars, research on the dynamic interaction between vehicles and road systems [3,4], and the vehicle–road dynamics model of coupled systems [5–11], our understanding in this area is continuously developing. The key problems that hinder the long-term development of heavy haul transportation have been partially solved, including the impact coefficient of vehicles, the vibration of vehicles under the uneven road surface, and the instantaneous response of road vibrations.

As mentioned above, the dynamic characteristics of the roadbed are regarded as one of the important components of the system. At present, a vehicle–road system can be divided into a vehicle vibration model and a road subsystem model. The vehicle is the moving part, which is the main excitation part that causes the vehicle and road vibration. The pavement and the following parts are the road subsystem. Due to the long-term bearing of the dynamic load of the upper vehicle, including the repeated load of heavy vehicle movement, subgrade diseases form easily. In particular, the deformation of flexible asphalt pavement is often affected by many factors such as traffic congestion, joints, load capacity, and repeated vehicle loads [12,13]. All these factors affect the normal service of the road and make the road surface irregular, seriously affecting the normal operation of the road and

**Citation:** Liang, B.; Xiao, J.; Shi, S. Establishment of an Eleven-Freedom-Degree Coupling Dynamic Model of Heavy Vehicle-Pavement. *Symmetry* **2022**, *14*, 250. https://doi.org/ 10.3390/sym14020250

Academic Editors: Yang Yang, Ying Lei, Xiaolin Meng, Jun Li and Juan Luis García Guirao

Received: 31 October 2021 Accepted: 11 January 2022 Published: 27 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

even damaging it. The above analysis shows that the damage of roads not only involves damage to the road surface but also the excessive deformation of the roadbed. Therefore, when analyzing the dynamic interaction between vehicle and road, the vehicle vibration excitation is the key problem, which is, therefore, one of the main concerns of this work.

Scholars have studied the dynamic problems of vehicle and highway subgrades from different perspectives [14,15]. However, few of the existing models consider the vehicle–road system as a whole, and most models analyze the road dynamic performance under the dynamic load of vehicles. The literature [16] describes the lateral dynamics model of vehicles under road environment excitation. The authors of [17,18] studied the dynamic response of the vehicle–pavement coupling system based on the nonlinear Timoshenko beam method and the vehicle–road dynamic response of the multi-degree-of-freedom vehicle with a double-layer, rectangular, thin plate. Zhang et al. [19] studied the road deformation and crack propagation path under the action of a quarter of vehicles using a modified two-parameter foundation plate. The above analysis fails to fully reflect the interaction of a vehicle–subgrade system in the process of vehicle moving, and this is another main objective of the present study.

In the present study, we establish a seven-degree-of-freedom vehicle dynamic model and use the harmonic superposition method to simulate the process of road random vibration. The dynamic interaction between vehicle and road is analyzed theoretically. Then, the validity of the model is verified by numerical analysis, and the dynamic response is analyzed. The organization of the paper is as follows. In Section 2, the dynamic model and equilibrium equation of heavy load vehicles are established. In Section 3, the coupling dynamic equation of heavy load vehicle–road is established. Then, the vehicle–road system is analyzed from multiple angles through numerical simulation, and the validity of the established equation is verified in Section 4.

### **2. Vehicle Dynamic Model and Equilibrium Equation**


*2.9. IoT Enabled System Architecture*

This system is regarded as a multiple-user access system, allowing numerous users to connect to the cloud simultaneously, as shown in Figure 3. There is only a single universal receiver shared by all users. An IoT system with cloud administration was created to classify obesity. Because it is a distributed system, the cloud is the best solution for a healthcare system that enables doctors to obtain data more easily. Our suggested IoT system comprises four key phases: (1) data collecting, (2) textual data classification, (3) diagnosis, and (4) user interface. Its goal is to lower disease rates through early detection of obesity.

**Figure 3.** Proposed IoT system architecture.

This figure demonstrates that the user is the origin of the entire process. With the web application interface, users engaging with the server and application interfaces are directly coupled. Therefore, when a user interacts with the web interface, a specific request is sent to the server. Upon receiving a request, the server examines it to determine what is the need of the user (obesity prediction, check his/her history, download report, or doctor's advice). Then the server will decide where to transmit the user's request after considering the needs of the user. Therefore, the server looks for an expert system that can handle the user's request and deliver the results. The server assigns the user's duties after identifying the optimal expert system. The user's task is inputted into the expert system as a string because the entire model is reliant on textual information, which is utilized to identify obesity in its early stages. After receiving a string input, the algorithm eliminates any extraneous words that are found during the prediction stage. After eliminating superfluous words, the user-provided data are used by the prediction engine to make predictions. Following the calculation of the outcome, the results are sent to the expert system. The server receives the results that the expert system collects. After obtaining the expert system's results, the server sends it to the web interface, where the user can access his/her results and move forward in light of the report.