Robust inverse parameter fitting of thermal properties from the laser-based Ångstrom method in the presence of measurement noise using physics-informed neural networks (PINNs)

$α$

Thermal diffusivity (⁠ $m 2 s − 1$⁠)

$ω$

Angular frequency of heating (⁠ $rad s − 1$⁠)

$ϕ$

Phase delay (⁠ $rad$⁠)

$ρ$

Density (kg m $− 3$⁠)

$A$

Amplitude of steady-periodic temperature oscillations (K)

$c p$

Specific heat (⁠ $J kg − 1 K − 1$⁠)

$f$

Frequency of heating (Hz)

$H$

Sample thickness (m)

$h$

Convective heat transfer coefficient (W m $− 2$ K $− 1$⁠)

$k x$

Thermal conductivity in the in-plane $x$-direction (W m $− 1$ K $− 1$⁠)

$k y$

Thermal conductivity in the in-plane $y$-direction (W m $− 1$ K $− 1$⁠)

$k z$

Thermal conductivity in the cross-plane $z$-direction (W m $− 1$ K $− 1$⁠)

$L$

Loss function

$n$

Number of grid points or pixels

$P$

Real part of complex temperature amplitude (K)

$q ″$

Heat flux (W m $− 2$⁠)

$Q$

Imaginary part of complex temperature amplitude (K)

$r$

Radius (m)

$SNR$

Signal-to-noise ratio

$T$

Temperature (K)

$t$

Time (s)

$w$

Loss weight

Heat generated in electronic devices, such as semiconductor chips and packages, must be dissipated to ensure reliability and operation below temperature limits. Heat must flow from within the package to the surrounding ambient or coolant through a series of thermal management components that may include thermal interface materials, heat spreaders, and heat sinks. Optimum thermal management is typically achieved by using materials with high thermal conductivity, k, so as to minimize the resistance to heat flow. Traditionally, high-thermal-conductivity metals with k values ranging from $∼100$–400 W m $− 1$ K $− 1$ have been used as heat spreaders. But with recent advances in the materials technology, there is a growing interest in engineered materials and composites, which can offer higher thermal conductivities than conventional metals.^1–3 However, such materials can often exhibit anisotropic thermal properties due to their composite nature or manufacturing process. As an example, naturally occurring isotropic graphite has a k of $∼50$ W m $− 1$ K $− 1$⁠, but synthetic graphite can have an in-plane k of up to 2000 W m $− 1$ K $− 1$⁠, while only $∼10$ W m $− 1$ K $− 1$ in the through-thickness direction.⁴ 

Although there is no standard technique to measure the in-plane k of anisoptropic materials, there are a couple of conventional techniques, such as the Ångstrom method⁵ and the laser-flash method,⁶ which have been adapted for such measurements. Both of these techniques are generally applicable for measuring k of “bulk” materials given an assumption of isotropic properties; to characterize the properties of anisotropic materials, these techniques rely on the fabrication of multiple samples cut along different orthotropic directions to measure the properties in these directions. Recent technique developments have sought thermal characterization approaches that can accurately measure in-plane anisotropic properties of materials in one measurement with a single sample. To bridge this gap, we recently introduced a new method for the measurement of in-plane thermal properties of isotropic and anisotropic films or sheets.⁷ 

Briefly, our method is based on the traditional Ångstrom’s method⁵ (for characterization of thermal diffusivity along one direction in thin and long rod-like materials), but extended for characterization in two dimensions to measure the in-plane thermal properties of films and sheets. This two-dimensional technique measures the steady-periodic temperature response of a material to periodic heating using a laser. Spatiotemporal temperatures are sensed using high-resolution infrared (IR) imaging. The two-dimensional discretized heat diffusion equation (along the two in-plane directions) in the frequency domain is evaluated throughout the spatial domain to extract the thermal properties of the material, using the steady-periodic response of the material. The amplitude of oscillations and the phase delay at each pixel location in the domain are calculated using a Fourier transform. This information is then used to perform an inverse parameter fitting to extract the thermal properties of the material using a numerical least-squares fitting algorithm. However, the requirement to numerically evaluate derivatives from the measurement data in this fitting approach makes the thermal property extraction potentially sensitive to noise in the measurement, especially when the magnitude of temperatures oscillations and phase delay are relatively low. In this work, we present an alternative method for inverse parameter fitting based on physics-informed neural networks (PINNs) to increase the tolerance to noise.

Physics-informed neural networks are deep learning frameworks that combine the efficiency of machine learning algorithms and the fundamental physics principles to solve a partial differential equation (PDE), or a system of PDEs, describing a physical system. They have gained immense popularity over conventional numerical modeling tools for performing certain tasks that benefit from the automatic differentiation capabilities of neural networks to evaluate higher order derivatives,⁸ thus avoiding the discretization errors typically encountered in numerical schemes. PINNs can be used to solve both forward and inverse problems. In the context of the current work, solving a forward problem could mean estimating the temperature fields for known or unknown boundary conditions, given temperature measurement data available at limited collocation points. On the other hand, this same trained model could be used to solve an inverse problem of deducing the material properties. For instance, Cai et al.⁹ demonstrate this twofold capability of PINNs in the case of Stefan phase-change problems, wherein limited temperature data measurements within the system are used within a PINNs framework to resolve the temperature distributions in each phase, while also inferring their respective thermal diffusivity values.

All measurement data is subject to some level of noise that arises from inaccuracy of sensors, measurement errors, and inherent variability of complex systems. PINNs, in their basic form or modified versions,^10–12 are particularly useful for making predictions/estimations that are robust despite high noise in the measurement data used for model training or evaluation. Garcia Inda et al.¹³ use PINNs to successfully reconstruct electrical properties from MREPT (magnetic resonance electrical property tomography) of noise-contaminated images. Oommen and Srinivasan¹⁴ demonstrate the utility of a basic PINNs model to solve inverse heat transfer problems in the case of rectangular pin fins (with different material properties and subjected to different boundary conditions). They show that a PINNs approach is faster and more robust to noise in the data, even compared to conventional machine learning techniques.

In this work, we employ PINNs to analyze data obtained from the 2D laser-based Ångstrom method. By doing so, we circumvent the challenges linked to numerical differentiation, and thereby demonstrate the robustness of PINNs approach for accurate inverse parameter estimation in the presence of noise that is added to the data. Our work thus presents an alternative data processing technique that complements the 2D laser-based Ångstrom method, enhancing its capability in measuring the thermal conductivity of anisotropic film-like materials. Section II of this paper overviews the 2D laser-based Ångstrom method, the experimental setup, and the least-squares method for inverse estimation of the in-plane thermal properties. Section III describes the methodology used to generate datasets for the purpose of this study and the PINNs-based approach for the inverse fitting of the thermal properties. Section IV compares the performance of the least-squares fitting approach and the PINNs approach under different levels of noise.

The 2D laser-based Ångstrom method, a metrology technique developed by the authors⁷ for the characterization of isotropic and anisotropic films and sheets, extends the principles of the traditional Ångstrom method of thermal diffusivity measurement, which was designed for thin and long rod-like materials (1D conduction), to two-dimensional conduction in thin films and sheets. The experimental setup is schematically represented in Fig. 1(a). A non-contact and stationary heat source, such as a focused laser beam or an IR-based LED light source, is incident at the center of the back side of the specimen, which is suspended over a metallic heat sink with a circular opening. For specimens that may be transparent to the wavelength of irradiation from the heat source, a thin and thermally black metallic circular disk may be attached to the underside of the sample to act as an absorber. The time-periodic heat source causes a periodic temperature oscillation in the specimen, $T(x,y,t)$⁠, which is measured from the top side using an IR camera after the specimen reaches a thermally steady-periodic state. Figures 1(b) and 1(d) show simulated (see Sec. III A) representative amplitudes of oscillation and associated phase delay in the suspended region of a hypothetical material with thermal conductivity of $k x$ = $k y$ = 10 W m $− 1$ K $− 1$⁠, a thickness of 500  $μ m$⁠, and at a heating frequency of 10 mHz. In real experiments, the data has noise and 1(c) and 1(e) illustrate the case data with noise added to the frequency domain data (see Sec. III B).

It must be noted that this technique assumes that the temperature gradients across the thickness of the specimen are negligible, relative to the gradients in the in-plane direction. This assumption can be realized by setting the frequency of heating such that the thermal penetration depth into the specimen exceeds the material thickness. This measurement is based on the principles of using the amplitude and phase lag of the steady-periodic temperature response of the material, and hence, the thermal property extraction is independent of the heat input or the boundary conditions at the periphery of specimen.

In the least-squares fitting technique, the second-order spatial partial derivatives of $P$ and $Q$ are calculated numerically at the spatial locations of each data point. Because experimental measurements involve the use of IR imaging, noise in the data has the potential to impact the extracted properties. The conversion of the data to the frequency domain using Fourier transforms eliminates much of the noise in the time series. However, when the second-order partial derivatives are calculated numerically, the effect of spatial noise is amplified. This is exacerbated for measurements with high spatial resolution, where the pixel-to-pixel distance is small, and for highly conductive materials, where the overall temperature rise for a given heat input is low. To reduce the effects of spatial noise, the spatial maps of the real and imaginary parts of the complex amplitude, $P$ and $Q$⁠, are smoothed by applying a square-shaped spatial convolution filter (filter2 in MATLAB), typically with a kernel size ranging from $5×5$ to $11×11$ pixels.

In practice, the transient temperature data to be used for inverse extraction of the thermal properties would be collected from an experiment. Herein, for purposes of assessing the newly developed fitting techniques, we instead generate transient temperature data from a numerically simulated experiment in COMSOL^TM Multiphysics. The simulated experiment is a numerical model that replicates the experimental system, consisting of a model geometry of the heat sink, specimen, and other associated boundary conditions. These simulated experiments have been extensively described and validated in the authors’ previous work.⁷ Using numerical experiments allows us to assess the inverse fitting methods against ground truth data, which can be generated for any hypothetical material properties and with varying levels of added signal noise.

The simulated experiment setup is shown in Fig. 2. A periodically varying temperature boundary condition is applied at the central $∼3$ mm diameter of the specimen, in the form of $T(t)= T ∞+ T a m p , max(1+sin(2πft))$⁠, where $T ∞$ is the ambient temperature, $T a m p , max$ is the maximum amplitude of oscillations, $f$ is the periodic heating frequency, and $t$ is time. The heat sink is assigned a fixed temperature boundary condition, $T s i n k$⁠.

In the simulated experiment, the properties of the hypothetical material are assumed to be known and include the thermal conductivity (⁠ $k x$⁠, $k y$⁠, $k z$⁠), density (⁠ $ρ$⁠), and specific heat (⁠ $c p$⁠). The simulation generates the transient thermal response of the material, and the recorded output is the transient surface temperature map, $T s i m(x,y,t)$⁠. These data are then used for extracting the thermal conductivity of the specimen, using either the least square fitting approach described in Sec. II (also described and validated in Ref. 7) or the PINN-based approach to be described in Sec. III C, without assuming any prior knowledge of the input $k x$ and $k y$ that were used to generate the temperature response.

This work aims to assess the advantage of PINN-based analysis under experimental conditions for which the measured temperature data are subject to a high level of noise. The data generated from the numerical simulations are noise-free. Therefore, noise is explicitly added to these synthetic data. Two types of noise are considered:

1. Noise addition in the time domain: White Gaussian noise is introduced to the transient temperature data (T) in the time domain (before taking Fourier transform). This case of noise addition simulates the impact of inherent noise in the infrared detector response over time.

2. Noise addition in the frequency domain: White Gaussian noise is added to the spatial maps of the complex amplitude components (P and Q) of the temperature signal after taking the Fourier transform. This case of noise addition simulates the impact of spatial variations in the emissivity of sample or pixel-to-pixel variations in the detector response.

The addition of Gaussian noise to the matrices is performed in MATLAB using the awgn function. The level of noise is quantified via the signal-to-noise ratio (SNR), with noise increasing as the SNR decreases. Noise is incorporated at SNR values of 1, 5, 10, 20, 30, 40, and 50, where SNR = 50 is the lowest noise level while SNR = 1 is the highest noise level considered in this study.

The overall objective of using the PINN framework for inverse parameter fitting is to minimize the loss function (⁠ $L$⁠), which takes into account both the training errors from the neural network and the underlying physics governing equations. Figure 3 shows the fully connected feed-forward neural network architecture that is employed. The network comprises an input layer with two neurons representing input variables, the $x$ and $y$ spatial coordinates, and an output layer with two neurons corresponding to $P$ and $Q$⁠. The outputs are then used to compute the residuals of the physics PDEs [Eqs. (3a) and (3b)] using automatic differentiation and solve the inverse problem of predicting $k x$⁠, $k y$⁠, and $h$⁠. The model is implemented in DeepXDE, a Python library¹⁶ and Google Colaboratory(https://colab.google/) services are used to execute the code for training the PINN model and extracting the inverse parameters.

The process of choosing model hyperparameters, such as depth (number of hidden layers), width (number of neurons in each hidden layer), learning rate, and loss weights, was informed by a preliminary set of training experiments. These initial trials involved systematically testing various hyperparameters one by one using numerically generated data with known inverse parameters (⁠ $k x$⁠, $k y$⁠, and $h$⁠) described in Sec. III A. The primary goal of these training experiments was to identify the optimal value for each hyperparameter under consideration (while keeping the others fixed) that resulted in the most accurate inverse parameter predictions and minimized the time needed for convergence to those values. Our approach involved optimizing each hyperparameter sequentially, using the optimized values from one set of experiments for the next set aimed at optimizing another hyperparameter.

After hyperparameter tuning, the final neural network architecture has a depth of 16 layers with each layer having a width of 20 neurons. The hidden layers use tanh activation functions and a Glorot Normal initializer. The network is trained using the Adam optimization algorithm with a batch size of 512. The activation function, the initializer, and the optimization algorithm have been chosen based on the typical approach employed in the literature for similar problems. The initial learning rate is set as $10 − 3$ and is then scheduled to decrease as the training progresses using an inverse decay algorithm with a decay rate of 0.2 every 1000 iterations.

In Eqs. (6)–(9), $N$ is the number of training data points (or the collocation points) obtained from the numerical simulations. This number is typically in the range of 5000–10 000 data points. Each term of the loss function has an associated weight (⁠ $w 1$⁠, $w 2$⁠, $w 3$⁠, and $w 4$⁠). Based on the hyperparameter tuning trials, these weights were assigned as 100, 100, 1, and 1, respectively. It is important for the loss terms to be balanced in the loss function. The volumetric heat capacity $ρ c p$ is typically of the order of $10 6$ J K $− 1$ m $− 3$⁠. We, therefore, normalize the input variables, $x$ and $y$⁠, by dividing them by the absorber disk diameter (⁠ $3× 10 − 3$ m). The training typically requires $∼$50 000–200 000 iterations before convergence of the fitted parameters is achieved.

To first demonstrate the effect of noise on the dataset to be input into the inverse parameter fitting algorithm, Fig. 4 shows the magnitude of the complex amplitude (⁠ $P 2+ Q 2$⁠) and phase delay [ $tan − 1( Q / P)$] of the temperature signal for an isotropic material of $k x$ = $k y$ = 10 W m $− 1$ K $− 1$⁠. A clean signal (no noise; solid red line) is plotted in comparison with several different noise-added signals. When noise (SNR = 1) is added to temperature data in the time domain, it is filtered to some extent by Fourier transform, which depends on the number of time-periodic cycles used. As shown in Fig. 4(a), the amplitude remains relatively unaffected by the time-domain added noise. However, noticeable deviation from the clean signal is apparent in the noisy data for the phase delay signal. The deviations are higher when using a lesser number of 5 measurement cycles (solid green circles) and increase at distances moving away from the center. The effect of the added noise diminishes with an increase in the number of 100 measurement cycles (solid blue circles), with these data nearly overlapping the clean signal. Similarly, when noise (SNR = 1) is added directly to P and Q (solid black circles), a significant difference in the phase delay is observed in Fig. 4(b). The amplitude in this case is also more affected by the noise compared to the case where the noise is added to the time domain. The difference is also apparent in the 2D maps of amplitude and phase delay as shown in Figs. 1(c) and 1(e), respectively.

The performance of the least-squares approach (without using any spatial convolution filter to smooth the data) vs the PINN approach for inverse fitting of the thermal conductivity is assessed under varying degrees of noise. Figures 5 and 6 show this comparison for noise added in the time domain for two isotropic samples [Fig. 5(a): $k x$ = $k y$ = 0.28 and (b) 10 W m $− 1$ K $− 1$] and one anisotropic sample [Fig. 6(a): $k x$ = 2 W m $− 1$ K $− 1$⁠, (b) $k y$ = 10 W m $− 1$ K $− 1$]. Notably, the PINN approach (black solid stars) uses data with just five time-periodic cycles and yields a very accurate prediction of thermal conductivity values for both isotropic and anisotropic samples across all noise levels, with errors in estimation remaining below 1%. In contrast, the least-squares fitting method (solid light blue, red, or green circles) yields higher error and is sensitive to the number of cycles used. When fitting five cycles of data (matching the number of cycles for the PINNs approach), the least-squares method performs very poorly, with > $∼85$% estimation error for SNR < $∼10$⁠. In general, using the least-squares approach, it is observed that the noise added in the time domain impacts the thermal conductivity estimations significantly when SNR < 30 and is, therefore, very sensitive to even small amounts of noise. This fitting method encounters challenges due to discretization errors stemming from utilizing numerical differentiation to compute second-order derivatives. In contrast, the PINN approach capitalizes on automatic differentiation, enabling the evaluation of derivatives via chain rule (backpropagation) while training the neural network, making this inverse fitting method robust even under high signal noise (up to SNR = 1 in this study).

The poor performance of the least-squares fitting method when noise is added to the time-domain data can be partially mitigated by increasing the number of temperature measurement cycles used in the analysis. The Fourier transform filters time-domain noise more effectively as the number of cycles increases. For instance, at SNR = 10, the estimation error is reduced from $∼83$% at 5 cycles to $∼17.5$% for 100 cycles. However, the approach is nevertheless susceptible to high noise levels and the least-squares fitting approach results in an estimation error of $∼50$%–70% in the worst-case scenario (SNR = 1). This demonstrates a twofold advantage of PINNs over the least-squares fitting approach: first, PINNs enable accurate extraction of the thermal conductivity from measurement data having very high levels of noise that would not otherwise be possible; second, even in cases of low or moderate noise, the PINN fitting approach offers a reduction in the number of cycles that must be measured to achieve an accurate inverse parameter extraction, thereby reducing the measurement time.

When introducing different levels of noise directly to P and Q in the frequency domain, the extracted thermal conductivity for the different fitting approaches are plotted in Fig. 7 for an isotropic thermal conductivity of $k x$ = $k y$ = 0.28 W/(m K). It is evident that the least-squares fitting approach (solid blue circles) generally fails to accurately extract the thermal conductivity when SNR < 50, whereas the PINN approach (hollow black stars) effectively fits the data for SNR > 20. Notably, PINN predictions start to falter for SNR < 20. This limitation is attributed to higher loss weights (100,100) assigned to the neural network losses compared to those for residual-based losses (1,1). These unequal loss weights lead to data overfitting and, consequently, the overfitting of noise. Adjusting the loss weights such that they are equal for both the neural network training losses and the PDE residual losses (1,1,1,1) helps in parameter fitting under extreme-noise situations. Implementing these equal loss weights (solid gray stars), the PINN approach is shown to accurately extract the thermal conductivity across all the noise levels. However, a drawback in using equal loss weights is that the PINN fitting approach requires more iterations to converge to the optimal values and, hence, higher computational time.

Taking the PINN approach with the original loss weights (100,100,1,1), we compare the fitting performance against the least-squares fitting approach for materials spanning different thermal conductivity values (0.1, 1, and 10 W m $− 1$ K $− 1$⁠) and anisotropy ratios (1, 10, and 100). The frequency and thickness of the specimen are fixed at 100 mHz and 100  $μ$m, respectively. For comparison, the fitting performance is evaluated for both clean data and noisy data (SNR = 10, added to the frequency domain). The left panel in Fig. 8 shows the thermal conductivity estimations by the two approaches when analyzing the data without noise. Both the methods work well and predict the thermal conductivity to within 5% error (shaded green) in most cases. In the cases, where one or both of $k x$ and $k y$ are 0.1 W m $− 1$ K $− 1$⁠, the estimation error increases to 35% and 30% for the least-squares and PINN approaches, respectively (corresponding to the yellow shading). This discrepancy is not associated with the fitting approach, but rather a physical limitation of the measurement under these conditions. Namely, this is attributed to the frequency of laser heating and sample thickness, which in this scenario, there is an insufficiently small thermal penetration depth for the lowest thermal conductivity to satisfy the assumption of 2D in-plane heat spreading in the specimen. As described in Ref. 7, the measurement heating frequency needs to be appropriately selected based on the extracted property measurement to ensure this penetration depth condition is satisfied.

The difference in the performance of the two fitting approaches is apparent when significant noise is added to data, shown for SNR = 10 in the right panel in Fig. 8. The least-squares fitting approach has extreme $>90$% estimation errors (shaded red). In contrast, the PINN approach fits the noisy data well and predicts the thermal properties with low errors similar to the clean data, excluding the cases of low thermal conductivity for the reasons already described. Overall, the PINN method performs much better compared to the least squared fitting method in the case of noisy data.

Although the PINN approach exhibits superior performance, it is not without limitations. Similar to many optimization routines, it was observed that the efficiency of the PINN approach depends on the initial guess values of the parameters being predicted; if these values deviate significantly from the ground truth values, it may require a large number of iterations and, consequently, substantial computational time for the inverse predictions to converge to specific values. This challenge is particularly evident in scenarios with high thermal conductivity (above 100 W m $− 1$ K $− 1$⁠).

The 2D laser-based Ångstrom method relies on inverse parameter fitting routines to extract the unknown thermal properties from spatiotemporal measurement data. A least squares fitting approach used in our previous work that developed this measurement technique relies on numerical differentiation to obtain second-order derivatives required for fitting to the governing heat diffusion equation. This approach is, therefore, susceptible to discretization errors, especially in the case of noisy data. In this work, we introduce an alternative inverse parameter fitting approach using physics-informed neural networks (PINNs).

The robustness of the PINN approach to recover the correct thermal properties is assessed through the introduction of time-domain or frequency-domain noise into numerically generated data. The PINN approach is robust down to very low signal-to-noise ratios (SNR = 1) with time-domain noise, with errors less than 1%. This greatly surpasses the performance of the least-squares fitting method, which fails to predict the thermal properties when the signal-to-noise ratio falls below 30. Furthermore, this robustness of the PINN approach is achieved by processing data using only 5 time-periodic measurement cycles vs the lower accuracy using 100 cycles for the least-squares fitting approach. This comparative performance of the inverse fitting approaches holds true for diverse specimen types encompassing both isotropic (low and high thermal conductivity) and anisotropic materials.

However, when noise is directly introduced in the frequency domain, discernible limits on the levels of noise tolerable by the PINN approach are observed, while the least-squares fitting method proves inadequate across all noise levels. The PINN method begins to have compromised prediction accuracy for SNR = 10 and below. This is attributed to the selected model architecture hyperparameters, namely, overfitting of the noise due to the unbalanced loss weights that give priority to fitting of the neural network outputs vs the physical governing equations. Selection of these parameters is a known challenge of such machine learning methods, and we, therefore, demonstrate the potential for improving the fitting accuracy by equalizing these loss weights under extreme-noise situations. This work highlights the potential of PINNs for extending the capability of this technique for characterizing a broader range of materials with higher accuracy owing to robustness in the inverse parameter fitting under practical levels of measurement noise.

Financial support for this work provided in part by members of the Cooling Technologies Research Center, a graduated National Science Foundation Industry/University Cooperative Research Center at Purdue University, is gratefully acknowledged. S.S. appreciates the financial support from the Adelberg Fellowship (awarded by the School of Mechanical Engineering at Purdue University). The authors would like to thank Ritwik V. Kulkarni, Rohan M. Dekate, and Pranay P. Nagrani, graduate researchers at Purdue University for their discussions and assistance with implementation of the PINN approach.

The authors have no conflicts to disclose.

Shanmukhi Sripada: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Aalok U. Gaitonde: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Justin A. Weibel: Conceptualization (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Amy M. Marconnet: Conceptualization (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.