Periodic oscillation for a class of in-host MERS-CoV infection model with CTL immune response

1.   Introduction

In medical history, the palpation approach is used to evaluate the pathological structure's stiffness to detect and classify pathology. In ancient Egypt, humans learned that a stiff structure in an organ (soft tissue) is an abnormal sign. Accordingly, palpation is used to detect and diagnose sickness; surgeons also applied this technique during operations and surgical intervention to detect diseased tissue. It can be said that the physical properties of tissues in terms of elasticity and viscosity provide helpful information in medical diagnostics, especially tumor detection ^[1]. According to ^[1], different types of tissue in a region (such as normal and diseased tissue (cirrhosis) of the liver tissue) give different values of elasticity. It is the basis for tissue elasticity estimation methods to diagnose the disease. Elastography includes two main steps: Firstly, it creates an excitation into the medium; secondly, it measures mechanical response from the medium and computes elasticity by different techniques. The stimulating method for shear wave propagation in tissues is classified by the stimulating method using acoustic radiation force or mechanical force. Elastic imaging techniques can also be classified according to point measurement methods or imaging methods. There have been many different techniques for creating elastic 1D, 2D, and 3D images. Shear wave elasticity imaging is introduced as a new ultrasonic technology of medical diagnostics in ^[2]. Complete principles and techniques of ultrasound elastography are then presented in ^[3]. In ^[4,5,6,7], the authors applied ultrasound elasticity imaging in soft mediums. In ^{[8,9,10,11,12]}, some advanced techniques were applied to improve the ability of ultrasound elastography.

Viscosity and elasticity estimations of soft tissues have created a significant leap in biomedical imaging and tissue classification. Shear wave elastography imaging (SWEI) was initially applied as extra medical diagnosis information ^[1]. Afterward, in 2013, Gennission et al. stated that SWEI provides significant advantages over other imaging techniques in quantitative imaging and elasticity contrast of tissues ^[13]. In 2014, Ferraioli et al. demonstrated applying elastography using shear waves to assess liver disease states ^[14]. Nowadays, SWEI equipment is commercialized ^[6]. However, equipment capable of imaging both elasticity and viscosity are rarely available.

To find these tissues ' characteristics, one needs to estimate shear wave propagation velocity in tissues at one or several frequencies. There are various ways to do so, such as acoustic radiation, surface vibration, or deep vibration as a mechanical force ^{[15,16,17,18,19]}. Orescanin et al. ^[20] obtained the shear wave propagation using a vibrating needle. Helmholtz algebra inversion transformation can be used to calculate the CSM directly ^[12]. However, this costs a lot of time and is easily affected by noise. In ^[21], the author analyzed errors of Helmholtz inversion for incompressible shear, vibration elastography. These drawbacks can be overcome by applying the Bayesian method with better velocity estimation as the low signal-to-noise ratio (SNR) ^[22]. However, this method is only applicable in a homogeneous medium of tissue. In our recent works ^[23,24], we have expanded the approach proposed in ^[22] into an inhomogeneous medium. Particularly, we applied the MLEF to efficiently estimate the CSM ^[25]. This method helps enhance the accuracy of CSM estimation while reducing the computational complexity. A limitation with this method is that we cannot detect the geometry properties of an object in the medium (i.e., a tumor in the body).

Most current methods are proposed for the automatic classification of soft tissues ^{[26,27,28,29]}. They are based on traditional techniques ^[30,31] and various machine learning algorithms such as decision trees ^[32,33], k-means ^[34], HMMs ^[35], and SVM ^[36]. These techniques have also been used to distinguish tissues, lesions, tumors, etc. Within the scope of this work, we emphasize the advantages of CSM estimation based on the extended Kalman filter (EKF) and Decision tree algorithm (DT) to categorize three types of soft tissue in the liver. Specifically, the proposed DT algorithm can differentiate between ordinary fibrosis, substantial fibrosis, and cirrhosis. Numerically simulation scenarios have been performed to demonstrate the effectiveness of the proposed EKF-DT method. The results illustrate that the EKF-DT approach improves the estimation quality and effectively classifies soft-tissue types.

2.   Methodology

2.1.   Detecting tumors in soft tissues

Figure 1 illustrates the imaging system for estimating shear wave velocity. The excitation source is a stainless needle vibrating at a specific frequency (100Hz ≤ f ≤ 500Hz). The shear wave afterward propagates in the tissue with the velocity of the particle measured by Doppler equipment.

In the mathematical view, the complex shear modulus (CSM) denotes quantitative information about the physical properties of soft tissues as below ^[1],

where ${\mu }_{1}$ is elasticity and ${\mu }_{2}$ is viscosity calculated by the Kelvin-Voigt model.

The region of interest (ROI) includes objects which need to be reconstructed (if present) in two-dimensional space (see Figure 2). The entire ROI is scanned over the angle from 0⁰ to 90⁰ with an angular shift of 1⁰ to create 91 rays. In each ray, the particle's velocity at each spatial point is collected by ultrasonic devices based on the Doppler effect. The following two objective functions are established to detect tumors by estimating the physical properties of tissues in terms of elasticity and viscosity:

and

For abnormal tissue, $\mu$₁ denotes elasticity, and ${\mu }_{2}$ denotes viscosity (if present); similarly, for ordinary tissue, ${\mu }_{1}^{0}$ and ${\mu }_{2}^{0}$ represents elasticity and viscosity, respectively, and $\omega$ denotes angular frequency ($\omega = 2\pi f)$.

In practice, it is impossible to estimate the CSM directly; instead, the CSM is computed through the wavenumber ${k}_{s}$ and attenuation coefficient α ^[13],

where $\rho$ is the density of the background medium ($\rho$ = 1000 kg/m³).

2.2.   Extended Kalman filter

The primary function of the Extended Kalman filter is to evaluate the wavenumber ${k}_{s}$ and attenuation coefficient α. At each spatial point, we develop the model by applying the extended Kalman problem.

Eq. (5) is used to calculate particle velocity $v(r, t)$ at position r, and time t ^[22]:

where A represents needle amplitude, ${r}_{0}$ denotes needle position, and $\phi$ indicates the initial phase; $\alpha$ and ${k}_{s}$, respectively denote attenuation coefficient and the wave number at spatial position r.

In discrete form, Eq. (5) is represented as below:

where $n$ is the time step-index and $∆t$ is the sampling period. Through the trigonometric transform of Eq. (6), we have:

To evaluate the wavenumber ${k}_{s}$ and attenuation coefficient α, Eq. (7) is described as state equation form:

where $f$ is a nonlinear vector function that describes the dynamics of the state vector.

Eq. (8) corresponds to

where ${x}_{n} = \left[\begin{array}{c}{v}_{n}\\ {\alpha }_{n}\\ {k}_{s\left(n\right)}\end{array}\right]$ is state vector at each point, the random variable ${p}_{n}$ is the process noise, ${k}_{s}$ is the wavenumber, ${k}_{s\left(n\right)}$ is the wave number evaluated at the time step n Because ${\alpha }_{n}$ and ${k}_{s\left(n\right)}$ are time independence; thus, we use $F$ as a nonlinear function that is the time dependence of ${v}_{n}$ at a certain position. $F$ describes the relation between ${v}_{n-1}$ and ${v}_{n}$ as shown in Eq. (7), ${\alpha }_{n-1} = {\alpha }_{n}$ and ${k}_{s(n-1)} = {k}_{s\left(n\right)}$ with supposing that α and ${k}_{s}$ are unchanged throughout the experiments. By using Doppler acquisition, the measured particle velocity at each spatial point is impacted by Gaussian noise ${\omega }_{n}\left(r\right)$. Thus, the measured particle velocity is:

To apply the EKF, Eq. (10) is rewritten in the measurement equation form of the Kalman problem:

Eq. (11) corresponds to

where ${{y}_{n} = \widehat{v}}_{n}$ is the measurement vector at each point. From Eqs. (8) and (11), ${x}_{n}$ is estimated by using the EKF corresponding to the algorithm in ^[30]. As a result, wave number ${k}_{s}$ and attenuation coefficient α at each point are estimated.

Afterward, we received a point cloud in ROI, where each point is presented by a pair of values, including estimated wave number ${k}_{s}$ and attenuation coefficient α. Therefore, the CSM coefficient is estimated by using Eq. (4). Finally, objective functions are reconstructed by using Eqs. (2) and (3).

For a quantitative evaluation of the effectiveness of 2D CSM estimation using the proposed method, a normalized error (or percent of estimated error, %RRE) was used and defined as follows:

where ${ϵ}_{\mu }$ and ${ϵ}_{\eta }$ are the normalized errors of elasticity and viscosity 2D estimations. $M\times N$ is the image size, ${\mu }_{i, j}$ and ${\widehat{\mu }}_{i, j}$ is the ideal elasticity, and estimated elasticity at pixel location $\left(i, j\right)$, ${\eta }_{i, j}$ and ${\widehat{\eta }}_{i, j}$ is the ideal viscosity and estimated viscosity at pixel location $\left(i, j\right)$. From that, we can calculate the ROC (receiver operating characteristic) value as 100%—%RRE. In the simulation scenario, Gaussian noise is generated to mix with the particle velocity of the shear wave. The degree of interference is expressed by the signal-to-noise ratio (SNR) parameter and is defined as follows:

where, P_signal is the signal's average power, and P_noise is the average power of the noise.

Figure 3 presents the proposed approach for two-dimensional CSM estimation and soft tissue classification using the Extended Kalman Filter (EKF) and Decision Tree (DT) algorithm, named the EKF-DT approach.

3.   Simulation and results

3.1.   The estimated results of the proposed EKF-DT algorithm

In this paper, simulation scenarios were implemented to verify the proposed EKF-DT method. The medium (soft tissue) has a size of 12.6 × 12.6 mm and three regions in it, corresponding to three pathological conditions of the liver. Three different tissue types with respective elasticity and viscosity are shown in Table 1. These values are referenced from work ^[26].

The CSM data was collected at 43 spatial locations per ray. The distance between two consecutive spatial positions was 0.3 mm, the sampling frequency was 10 kHz, 500 samples were collected at each spatial position, and the medium density $\rho$ = 1000 kg/m³. The vibration needle had a diameter of 1.5 mm, and the vibration frequency of the needle that we chose is 100Hz. For the frequency value f, essentially, a small value is good because it can penetrate tissue deep enough to image tissue far from the stimulus. However, if a too low frequency is chosen, the shear wave's frequency oscillation may be equal to or slightly larger than the disturbance frequency. Then, the estimated quality will be very poor. Therefore, the shear wave frequency should not be too low and should not be too high. In this paper, Threshold_A and Threshold_B are initially selected as (2.310 kPa, 1.885 Pa.s) and (3.620 kPa, 3.146 Pa.s), respectively. The steps to measure the particle velocity of the shear wave in the ray scanning method are as follows: First, the evaluated area of tissue is divided by 91 rays evenly spaced along the angle from 0⁰ to 90⁰ in polar coordinates, where the coordinate (0, 0) is also the position of the vibrating needle. Next, a shear wave generator at a frequency of 100 Hz is excited at this particular location through a vibrating needle with the amplitude A = 8.6 mm. Finally, the linear array transducer BW-14/60 was used to measure the particle velocity of the shear wave at 43 equally spaced locations on each ray. In total, we need to measure the particle velocity at 91 × 43 = 3913 locations.

Figure 4 illustrates the objective function of ideal elasticity O₁, where the vibration needle is located at the origin of coordinate O. The values of the three tissue types are shown in Table 1. The elasticity discrepancy between fibrosis and healthy liver is not noticeably represented. The ideal initialization image (Figure 4) has a size of 12.6 × 12.6 mm. The vibrating needle is located at the upper left vertex of the rectangle. Centered as the needle placement, the region of interest (ROI) will be an arc with a radius R = 12.6 mm. The ROI contains three regions corresponding to three pathological conditions of liver tissue (normal, fibrosis, cirrhosis).

Figure 5 illustrates the recovered elasticity objective function ${\widehat{O}}_{1}$ using the EKF-DT method. Despite the appearance of ripples on the surface, the reconstructed results closely resemble the ideal function. These ripples result from additive noise and the imperfection of the model. It is not easy to differentiate between fibrosis and cirrhosis because of ripple influence. Moreover, Figure 5 indicates that the ray scanning technique cannot sweep through the entire square area of 12.6×12.6 mm². The region of interest (ROI) is limited to a quarter circle area with a radius of 12.6 mm. Due to the hardware limitation of the BW-14/60 linear array transducer, it is possible to measure the particle velocities of 43 equally spaced locations per ray. Therefore, the scannable length is also the radius of the circle centered at the vibrating needle. Thus, points outside the scan radius cannot be estimated by the limitation of this transducer. In fact, to investigate a larger area, it is necessary to slide the transducer; the final image is a synthetized image of the obtained small images.

Figures 6a and 6b shows the estimated wave number and attenuation coefficient along the 60th ray with SNR = 30dB. It is noticeable that a casual change of the wavenumber happens at the distance of 6 mm (between fibrosis and cirrhosis) and 100 mm (between cirrhosis and healthy liver tissue) compared to the coordinate origin O. Similarly, this change can be seen in Figures 4 and 5. Figures 6a and 6b show that the estimated wave number and attenuation coefficient can fairly approach the ideal lines. Nevertheless, the ripples are significantly more significant than the estimated wave number when evaluating the attenuation coefficient. After assessing the wavenumber and the attenuation coefficient, the elasticity and viscosity are calculated using Eq. (4). Finally, the objective functions O₁(r) and O₂(r) are reconstructed.

Applying Threshold_A and Threshold_B in the decision tree algorithm (DT), as seen in Figure 3, can easily separate cirrhosis locations from the fibrosis and healthy tissue. Figure 5 indicates a reconstructed CSM image that accurately shows the location of soft tissues (fibrosis, cirrhosis, and healthy liver) in the ROI.

3.2.   Find the best values of thresholds

3.2.1.   Investigate the thresholds for elasticity

The selection of the initial threshold pair for elasticity is $\left({\mu }_{1T} = {\mu }_{11}+\frac{{\mu }_{12}-{\mu }_{11}}{2}, {\mu }_{1{T}'} = {\mu }_{12}+\frac{{\mu }_{13}-{\mu }_{12}}{2}\right)$ = (2310, 3620). To investigate the increase of this threshold pair to value pair (2560, 4680), we take ten investigated points, then the hopping value pair is (25,106). Similarly, to investigate the decrease of this threshold pair to the value pair (2060, 2560), we also take ten investigated points, and the hopping value pair is also (25,106).

Figure 7 shows the influence of different threshold pairs of elasticity and viscosity on the estimation quality, and Figure 8 shows the ROC (receiver operating characteristic) curve to obtain Threshold_A and Threshold_B. From Table 2 and Figures 7–8, it can be seen that the smallest estimated error percentage of elasticity is 2.6982% with the threshold pair of (2310, 3620).

3.2.2.   Investigate the thresholds for viscosity

The selection of the initial threshold pair for viscosity is$\left({\mu }_{2T} = {\mu }_{21}+\frac{{\mu }_{22}-{\mu }_{21}}{2}, {\mu }_{2{T}'} = {\mu }_{22}+\frac{{\mu }_{23}-{\mu }_{22}}{2}\right)$ = (1.995, 3.73). To investigate the increase of this threshold pair to value pair (2.27, 5.19), we take ten investigated points, then the hopping value pair is (0.0275, 0.146). Similarly, to investigate the decrease of this threshold pair to the value pair (1.72, 2.27), we also take ten investigated points, and the hopping value pair is also (0.0275, 0.146).

Table 3 and Figures 7–8 show that the smallest estimated error percentage of viscosity is 2.1834%, with the threshold pair of (1.885, 3.146).

After well-investigation, Threshold_A and Threshold_B are found as (2.310 kPa, 1.885 Pa.s) and (3.620 kPa, 3.146 Pa.s), respectively. These pairs of thresholds can be used to obtain the best estimation performance with the smallest estimated error percentage.

4.   Discussions

To estimate a CSM image of soft tissue elasticity and viscosity, we used an extended Kalman filter (EKF). Then, based on the estimated CSM image, we apply the decision tree (DT) algorithm to classify three pathological conditions of the liver (normal, substantial fibrosis, and cirrhosis).

First of all, on the adoption of an EKF filter, this usage has the following advantages. For the first superiority of the EKF scheme, we just use a dataset obtained at a single oscillation frequency (100 Hz) in the absence of fitting multiple-frequency information to a rheological model in the viscoelastic environment. Correspondingly, the acquisition time can be remarkably saved in comparison with the state-of-the-art approaches ^[37,38]. Thereto, this EKF approach can work with a high frequency as 100 Hz, meanwhile, the approaches in ^[37,38], being oscillatory or dynamic tests that may be examined as frequency domain solutions, only have the ability to study the viscoelastic characteristics at low frequencies (smaller than 100 Hz) because inertial effects cause biasing with measurements at higher frequencies ^[39]. For the second superiority, different from the former approaches in ^[40,41], the EKF scheme encompasses noise reduction due to utilizing the boon of the EKF. The fact that the noisy influence is a heavy challenge in realistic applications. For the third superiority, the EKF's low complication is appropriate for estimating the CSM in real-time or online. As analyzed in ^[42], the imaging total runtime of the EKF is smaller than the one of the MLEF. Notwithstanding, the proposed EKF scheme only concern about estimating the one-dimensional CSM and imaging of two-dimensional visco-elasticity of soft tissues, while it is a more practical issue of three-dimensional viscoelastic imaging. This approach deals with the problem of tracking two-dimensional shear-wave, also known as tracking in-plane wave propagation. And so, it cannot accurately track out-of-plane wave propagation, causing an erroneous CSM estimation. In addition to that, the approach works well for the simulation data, but experimental data is of value and is an issue for later consideration. Ultimately, the CSM estimation with the EKF scheme is contingent on the tissue response dynamic model, which is a constitutive viscoelastic model. Regularly used constitutive models are idealized, such as Maxwell, Voigt, and Kelvin models, and do not directly use for natural tissues ^[43]. It is still an open research issue.

To discuss algorithms for classifying tissue pathological conditions, four algorithms were examined in ^[44], with paying attention to expressing the traditional statistical spectrum and late progress in schemes based computer: discriminant analysis; regression models; tree-based algorithms; artificial neural networks. Each scheme is assessed according to four criteria: precision needed computing time for producing results, understandability of the results, and the easy use of the algorithm to relatively naive medical users. In this paper, we choose tree-based algorithms because this algorithm has no hypothesis regarding the underlying distribution; therefore, this approach executed consistently well. Indeed, in this paper, based on the CSM estimation results of the EKF algorithm, we applied the decision tree algorithm to classify well the three histological conditions of the liver with the performances of elasticity and viscosity classification are 97.3018% and 97.8166%, respectively. Moreover, with the best estimation performance, we find out the optimal classification threshold pair of (2.310 kPa, 1.885 Pa.s) and (3.620 kPa, 3.146 Pa.s), respectively. The number of locations will become much larger depending on the estimated image resolution and array probe technology. And we have to identify each point in the processing space, leading to a large number of estimated points, so a fast classification algorithm is needed to perform. In addition, the number of features in the problem is not much, so a simple algorithm like a decision tree is very effective.

Several other effective solutions have also been studied to classify the disease states of the tissue, such as ^[45,46,47]. In ^[45], a deep learning (DL) structure is built to automatically extract the features of the learned-from-data image in the shear wave elastography, and the DL architecture differentiated benign and malignant breast tumors. In ^[46], with integrating two quantitative ultrasound parameters, SWVmax and RI, a plain classification scheme represents a high diagnosis performance. This scheme can exactly discriminate between benign and malignant breast lesions. Furthermore, this can reduce the unessential number of breast biopsies up to 60 percentages for all cases, dodging the bias in subjective explanation. In ^[47], an approach combing of ultrasound, shear wave elastography, and MRI is proposed to predict and characterize benign or malignant tissue precisely. With the inclusion of shear wave elastography, it is supposed to specify the threshold value of velocity for classifying benign or malignant lesions.

With the three advantages of EKF filter: only need to collect data set at a specified vibration frequency without using multi-frequency information; exploit noise reduction to the advantage of EKF; low complexity of EKF is suitable for real-time/online CSM image imaging, this filter can be used for real-time/online image imaging. In ^[48], the authors used an adaptive filter that faced the real-time requirement. Currently, with the global Covid-19 epidemic situation, remote/online diagnosis of medical conditions is essential and urgent. Therefore, based on the advantages of the EKF filter, we continue to apply a simple machine learning algorithm (decision tree algorithm) to classify the pathological condition of tissues to support the remote diagnosis, supporting the doctor's diagnosis. However, the proposed approach also has some limitations as follows: (a) To be able to apply the EKF filter, the particle velocity of the shear wave is modeled according to equation (5) which is only suitable for homogeneous environments; (b) Imaging by ray scanning is quite simple because it can be converted to the CSM estimation at each tissue space location. The limitation of this method is that it only estimates CSM in the ray-scanning region, whose center is the excitation site, and the estimated locations are not equally spaced (the region far from the excitation source will be sparser than the region near the excitation source); (c) Decision tree algorithm is still quite simple, it needs to be developed or applied by more complex machine learning algorithms; (d) The study is just implemented with the simulation data set, to be able to apply in practice, the experimental implementation needs to be carried out.

5.   Conclusion

This paper has demonstrated that the CSM estimation of soft tissues could find and classify several significant states of tissues. Numerically simulated scenarios reconstructing objective functions of elasticity and viscosity have been implemented to confirm the EKF-DT method's efficient outcomes. The simplified DT algorithm efficaciously classifies the states of estimated tissues. As a consequence, this approach can automatically spot health issues by categorizing three states of the liver (normal, substantial fibrosis, and cirrhosis) with the well-investigated threshold pair. In the next work, we improve classification performance by integrating the DT algorithm with SVM. Moreover, 3D image recovery and experiments could be used to verify the suggested EKF-DT solution ^[49,50].

Funding

This work was supported by the Ministry of Education and Training under Grant No. B2020-SP2-02.

Conflict of interest

The author declares that he has no conflict of interest.