The results of several scientific and clinical studies on AF have led to the development of many mechanistic and nonmechanistic approaches for treating AF. Every treatment proce‐ dure has its own benefits and limitations. So, the treatment method that gives the maximum patient outcome is what is followed widely. Each treatment also has its specific endpoint, however, the ultimate goal of any treatment is to successfully terminate AF and restore sinus rhythm. Some common treatment methods are discussed below.

#### *2.3.1. Pharmocological therapy*

The studies on different ionic mechanisms of the cardiac cells and how they are related to the arrhythmia (e.g., formation of AP alternans), have led to the development and synthesis of various rate‐control and rhythm‐control drugs. They directly interact with the ion channels, and control and regulate the exchanges in order to maintain a normal AP thus maintain‐ ing NSR. The antiarrhythmic drugs are categorized into various classes based on the behav‐ ior of the drug on different ion channels. Class‐I drugs are related to sodium channels, e.g., quinidine. Class‐II are effective beta blockers, e.g., propranolol. Class‐III type interacts with potassium channels and prolongs the action potentials, e.g., sotalol, and Class‐IV helps in regulating the calcium channels, e.g., diltiazem. Class‐III is a commonly prescribed antiar‐ rhythmic for the fact that it prolongs the AP, increases the refractory period, and so the AP cannot be fired instantly, thus preventing reentries. Some popular Class‐III drugs are amioda‐ rone and dofetilide. Unfortunately, the drugs may not be much effective as a standalone AF treatment especially for patients showing poor response [38], and it might also produce some side effects for the patients.

## *2.3.2. Direct current electrical cardioversion*

Direct current electrical (DCE) cardioversion or simply cardioversion, is the procedure by which an electric voltage is applied to the chest at particular instants of the heart beat in order to reset the rhythm back to NSR [39]. The electric shock at a specific power is applied to the chest using electrode patches to the patient under anesthesia. The power is decided according to the intensity of AF. The shock is provided in such a way that it synchronizes with the R‐wave of the PQRS complex, mainly because the atria would have passed the refrac‐ tory period during R‐wave, and the application of the current will excite the cells and start a fresh cardiac cycle. This treatment, most of the time terminates AF immediately, however, some patients experience AF episodes return after cardioversion. Hence, it is important for the clinician to make a reliable decision on whether or not the patient would benefit from DCE cardioversion therapy.

#### *2.3.3. Open‐heart maze procedure*

The maze procedure is a surgical treatment by which lesions are created using incisions to various regions of the myocardium after opening up the chest. The treatment technique was introduced by James Cox in 1991 [40]. The lesion set resembles the structure of a maze, which gives it the name. The lesions basically are designed in such a way that the signals in the atrium do not form a reentry. The maze procedure has proved to be effective in terminating AF; however, it is highly invasive and is normally provided for patients who undergo an open‐heart surgery as part of their treatment for other types of heart complications.

#### *2.3.4. Catheter ablation therapy*

Catheter ablation therapy is a method by which ablation lesions are formed using radio fre‐ quency (RF) energy delivered via a catheter. The conventional catheter ablation procedure for AF today is the pulmonary vein isolation (PVI). In PVI, the lesions are formed around the PV ostia with the aim of electrically isolating the PVs from the rest of the left atrium. The PVI came into existence because of the pioneering discovery by Haissaguerre et al. that some ectopic trig‐ gers from the PVs enter the atrium and maintains AF [41] (illustrated in **Figure 5**). Some studies also suggest that these focal triggers could initiate non‐PV AF drivers as shown in **Figure 5**. More about these non‐PV drivers are discussed in Section 3. PVI showed remarkable increase in suc‐ cess rate of AF treatment and hence was adopted as the standard procedure by clinicians across the globe. There was also improvement over the traditional PVI that proposed only lesions around PV ostia. A strategy called the circumferential ablation was introduced which defined the lesions around the antrum of the PVs as opposed to the ostia [42]. In addition to circumfer‐ ential ablation, a linear lesion set was proposed, which appeared to produce better outcome in clinical studies. Some of these non‐PV lesions consisted of roof line, lines at mitral isthmus, and lines in between the PV pairs. PVI is mostly the primary strategy, and linear lesions are mostly preferred as a strategy in addition to PVI. However, what ablation to perform on top of PVI is disputable and is decided mostly based on the patient's individual atrial and AF properties.

Catheter ablation treatment combined with antiarrhythmic drug therapy is followed widely today. However, the treatment is not completely successful. Many patients experience recurrence of AF after the first procedure. Multiple repetitions of the procedure have to be performed especially for persistent AF patients, to control the recurrence. After a single pro‐

**Figure 5.** PV triggers initiating AF (posterior view): Figure illustrates the ectopic focal triggers originating from the PVs eventually initiating spiral waves in the LA wall, thus initiating and sustaining AF. LSPV: left superior pulmonary vein; RSPV: right superior pulmonary vein; LIPV: left inferior pulmonary vein; RIPV: right inferior pulmonary vein.

cedure, late recurrence of AF occurs between 11 and 29%, and after repeated procedures, it occurs between 7 and 24% [37]. A widely accepted reason for the recurrence is the reconnec‐ tion of the PVs, however, there are also theories suggesting that the recurrence could be due to the arrhythmogenic sources (especially functional sources) existing outside the PVs, which are not addressed by the current catheter ablation procedure.

## **3. Computer‐aided technologies for atrial fibrillation**

AF; however, it is highly invasive and is normally provided for patients who undergo an

Catheter ablation therapy is a method by which ablation lesions are formed using radio fre‐ quency (RF) energy delivered via a catheter. The conventional catheter ablation procedure for AF today is the pulmonary vein isolation (PVI). In PVI, the lesions are formed around the PV ostia with the aim of electrically isolating the PVs from the rest of the left atrium. The PVI came into existence because of the pioneering discovery by Haissaguerre et al. that some ectopic trig‐ gers from the PVs enter the atrium and maintains AF [41] (illustrated in **Figure 5**). Some studies also suggest that these focal triggers could initiate non‐PV AF drivers as shown in **Figure 5**. More about these non‐PV drivers are discussed in Section 3. PVI showed remarkable increase in suc‐ cess rate of AF treatment and hence was adopted as the standard procedure by clinicians across the globe. There was also improvement over the traditional PVI that proposed only lesions around PV ostia. A strategy called the circumferential ablation was introduced which defined the lesions around the antrum of the PVs as opposed to the ostia [42]. In addition to circumfer‐ ential ablation, a linear lesion set was proposed, which appeared to produce better outcome in clinical studies. Some of these non‐PV lesions consisted of roof line, lines at mitral isthmus, and lines in between the PV pairs. PVI is mostly the primary strategy, and linear lesions are mostly preferred as a strategy in addition to PVI. However, what ablation to perform on top of PVI is disputable and is decided mostly based on the patient's individual atrial and AF properties.

Catheter ablation treatment combined with antiarrhythmic drug therapy is followed widely today. However, the treatment is not completely successful. Many patients experience recurrence of AF after the first procedure. Multiple repetitions of the procedure have to be performed especially for persistent AF patients, to control the recurrence. After a single pro‐

**Figure 5.** PV triggers initiating AF (posterior view): Figure illustrates the ectopic focal triggers originating from the PVs eventually initiating spiral waves in the LA wall, thus initiating and sustaining AF. LSPV: left superior pulmonary vein; RSPV: right superior pulmonary vein; LIPV: left inferior pulmonary vein; RIPV: right inferior pulmonary vein.

open‐heart surgery as part of their treatment for other types of heart complications.

*2.3.4. Catheter ablation therapy*

100 Computer-aided Technologies - Applications in Engineering and Medicine

The clinical DSS is largely used during an AF ablation procedure. In order to ablate different regions of the LA, it is necessary to have a visualization of both the atrium and the catheter's location. The fluoroscopy technique which produces a real‐time image of the catheter and the atrial structures was used initially for visualization. However, fluoroscopy uses X‐rays, and prolonged exposure is harmful for the patient as well as for the clinicians performing the procedure. In order to reduce the fluoroscopy time, a novel technology called the electroana‐ tomic mapping was introduced. Electroanatomic mapping is a technique where the clinician maneuvers a mapping catheter across the LA endocardium and the computer produces a three‐dimensional virtual map of the LA endocardium. This mesh structure displays the LA anatomic structures along with the PVs (the atrial structure shown in **Figure 5** is a mesh pro‐ duced by such a mapping system). In addition, the electroanatomic mapping systems have the function to display the real‐time catheter position on the mesh. This makes it possible for the clinician to ablate anywhere using the 3D mesh as the reference. Some electroanatomic mapping systems available in the market are "CARTO" by Biosense Webster Inc., "EnSite NaVx" by St. Jude Medical, and "Rhythmia" by Boston Scientific. Electroanatomic mapping systems not only produce the mesh structure, but also are capable of performing several other interesting functions that help clinicians to make ablation decisions. One such feature is the activation map. As the name "electroanatomic mapping" suggests, the map shows anatomy as well as the electrical activations in the form of a voltage map with color codes. This largely helps in treating atrial tachycardia, the clinician can see in the electroanatomic map and visu‐ alize where the earliest activation occurred, and then perform ablation accordingly. Further, the system displays the ablation lesions on the mesh with a unique color code that is varied based on the power delivered by the RF catheter. This helps the clinician to ensure that the tissue was sufficiently cauterized to stop the conduction at that location.

The mapping catheters that are used to construct the electroanatomic mesh are special‐ ized diagnostic catheters known as a multipole diagnostic catheter (MPDC) (some types of MPDCs are shown in **Figure 6**). It consists of multiple nodes with voltage sensitive electrodes. These electrodes are arranged in a particular geometric shape which differs between differ‐ ent MPDC types. Some common MPDCs are Lasso and Pentaray both by Biosense Webster, FirMap (Abbott EP), Constellation (Boston scientific), Achieve (Medtronic), and Advisor (St. Jude Medical). A major feature of the MPDCs is that the electrodes can record the EGMs when they are placed in contact with the tissue. Hence, the MPDC is placed in different endocardial locations to build the mesh structure and also to obtain the EGMs. These EGMs can further be used to perform non‐PV ablation. As stated earlier, PVI is not completely successful, and in order to improve the success, non‐PV targets should be determined and ablated. For example,

**Figure 6.** Different types of MPDCs: Figure shows the various types of MPDCs, which are used for generating electroanatomic maps and for recording information from the endocardium. The top one is a 20 unipole Lasso catheter, the bottom‐left is a Basket catheter and the bottom‐right one is a Pentaray catheter. Image reproduced from Ref. [33].

some studies showed the presence of slow conduction regions in the myocardium, which when ablated, terminated AF in their studies [43]. These sites can be identified by recording the EGMs, which, based on a specific definition, could be characterized as what is called as complex fractionated atrial electrograms (CFAEs). CFAE ablation was the first EGM‐based (mechanistic) non‐PV ablation, while all other strategies are anatomy based (nonmechanistic), in the sense that they do not require any functional quantities such as the EGMs to determine whether or not to ablate that site. Based on the above concept, some mapping systems were installed with a feature to automatically compute CFAE sites, and indicate them on the mesh, for the clinicians to make the ablation decision about CFAE. However, in current clinical practice, CFAE ablation is not very common, possibly because of some controversial results reported by some studies [37, 44]. This demonstrates the need for determining the "real" non‐ PV targets which play a role in AF maintenance. Traditionally, in AF studies performed in the past using physical experiments and computer simulations, several studies have demon‐ strated circular and spiral patterns of electrical waves propagating along the tissue during AF [45–48]. The AF spiral waves were termed as rotors which later gained more attention among researchers in the field. Rotors consist of a wavefront and a wavetail, and they may rotate about a functional center (**Figure 7**). The gap between the front and the tail is the wavelength of the rotor. A rotor propagates with a conduction velocity, and rotates following a peri‐ odic cycle length (the velocity vectors at the wavefront of a rotor are illustrated in **Figure 7**). Rotors can also form wave breaks that eventually lead to multiple chaotic spiral waves. Rotors were initially encountered only in animal studies. For the first time, they were claimed to be observed in human clinical studies by Narayan et al. [49]. In addition to finding rotors, their study also reports observation of some ectopic foci outside the PVs, and that they act as AF drivers. The technique that revealed the presence of rotors using the EGM in humans is called phase mapping [49]. According to this technique, the EGMs are collected simultaneously at

**Figure 7.** Propagation of a rotor wave: The different components of a rotor wave are illustrated in the figure. The arrows points in the direction of the wavefront propagation. Image recreated based on the work in reference [50].

some studies showed the presence of slow conduction regions in the myocardium, which when ablated, terminated AF in their studies [43]. These sites can be identified by recording the EGMs, which, based on a specific definition, could be characterized as what is called as complex fractionated atrial electrograms (CFAEs). CFAE ablation was the first EGM‐based (mechanistic) non‐PV ablation, while all other strategies are anatomy based (nonmechanistic), in the sense that they do not require any functional quantities such as the EGMs to determine whether or not to ablate that site. Based on the above concept, some mapping systems were installed with a feature to automatically compute CFAE sites, and indicate them on the mesh, for the clinicians to make the ablation decision about CFAE. However, in current clinical practice, CFAE ablation is not very common, possibly because of some controversial results reported by some studies [37, 44]. This demonstrates the need for determining the "real" non‐ PV targets which play a role in AF maintenance. Traditionally, in AF studies performed in the past using physical experiments and computer simulations, several studies have demon‐ strated circular and spiral patterns of electrical waves propagating along the tissue during AF [45–48]. The AF spiral waves were termed as rotors which later gained more attention among researchers in the field. Rotors consist of a wavefront and a wavetail, and they may rotate about a functional center (**Figure 7**). The gap between the front and the tail is the wavelength of the rotor. A rotor propagates with a conduction velocity, and rotates following a peri‐ odic cycle length (the velocity vectors at the wavefront of a rotor are illustrated in **Figure 7**). Rotors can also form wave breaks that eventually lead to multiple chaotic spiral waves. Rotors were initially encountered only in animal studies. For the first time, they were claimed to be observed in human clinical studies by Narayan et al. [49]. In addition to finding rotors, their study also reports observation of some ectopic foci outside the PVs, and that they act as AF drivers. The technique that revealed the presence of rotors using the EGM in humans is called phase mapping [49]. According to this technique, the EGMs are collected simultaneously at

**Figure 6.** Different types of MPDCs: Figure shows the various types of MPDCs, which are used for generating electroanatomic maps and for recording information from the endocardium. The top one is a 20 unipole Lasso catheter, the bottom‐left is a Basket catheter and the bottom‐right one is a Pentaray catheter. Image reproduced from Ref. [33].

102 Computer-aided Technologies - Applications in Engineering and Medicine

different LA locations using a basket type MPDC (shown in **Figure 6**), and the signals are processed using Hilbert transform in order to obtain a map of the spatiotemporal variation of the signal phase. This phase map showed progressive change of phase in a circular pattern around a core called as the phase singular, thus demonstrating the characteristic of a rotor. The method of Hilbert transform that is used to obtain the phase map has also been used in different contexts of fibrillation in earlier studies [51, 52]. The algorithm/ablation procedure used by Narayan et al. was called focal impulse and rotor modulation (FIRM). Based on the outcome of the clinical studies on FIRM‐guided ablation and several other follow‐up studies, the FIRM algorithm was found to be effective, especially when combined with PVI [53]. The FIRM technology was later commercialized as a clinical diagnosis system for AF, now owned by Abbott Electrophysiology.

While whether or not the rotors are the primary AF drivers, is a topic still under investigation, developing novel algorithms for identifying and localizing rotors have been a critical area of research in the community. Several methods were introduced to help map rotors using engi‐ neering techniques in both time and frequency domain [54–58]. Some of these methods are specifically designed to use conventional MPDCs such as a Lasso catheter, aiming to outper‐ form the limitations of the current EP systems available for identifying rotors.

As opposed to the time domain approach of phase mapping, there were also some fre‐ quency domain approaches, one of which is called the dominant frequency (DF) mapping. According to the DF mapping technique, the fast Fourier transform (FFT) of the signals at various locations is computed, and the highest frequency from the FFT plot at each location is considered as the DF for the respective location. When the DFs at various regions of the atrium were investigated, the frequency distribution showed organization and was said to have high values in the left atrium [50]. However, there are some studies suggesting that DF is not an effective method to analyze AF [37].

The workstations available today in the EP labs consists of all the components required to assist the clinicians during an ablation procedure, such as the electroanatomic mapping soft‐ ware, mapping and ablation catheters, imaging devices such as fluoroscopy and ultrasound, and other associated software such as the FIRM mapping, DF mapping, and CFAE mapping. Therefore, EP workstation is the best example of a clinical DSS, without which the ablation procedures is nearly impossible to perform today. This demonstrates the significance of clini‐ cal DSS not only in AF, but also in general clinical cardiac electrophysiology.

## **3.1. Developed methods of clinical decision making for AF**

As discussed earlier, clinical DSSs play an important role in AF ablation treatment. This sec‐ tion presents a detailed discussion of some techniques and algorithms developed by our group that could be potentially implemented in a clinical DSS to assist a clinician for manag‐ ing AF patients. A schematic diagram of the various EP systems along with the integration of one of the developed clinical DSS system is illustrated in **Figure 8**.

The role of clinical DSS in AF could be categorized into three important tasks:

The above three tasks quite covers the necessary tasks required for AF management today. All the algorithms mentioned below were implemented in MATLAB (The MathWorks, Inc., Natick, MA, USA).

### *3.1.1. Detection of AF from surface ECG for AF therapy management*

Rate control and rhythm control are two strategies for treatment of AF. The goal of rhythm control is to convert an arrhythmic condition to NSR, and that of rate control is to slow down the heart rate when it increases due to AF. Catheter ablation therapy is an example of rhythm control method. Although catheter ablation is a commonly preferred treatment, at an initial stage of AF, doctors recommend rate control therapy. It could be performed using rate control drugs or cardiac pacemakers.

Detection of AF is a crucial operation for providing timely and appropriate treatment for the patient. AF detection techniques mostly involve application of signal processing algorithms to ECG signals [59]. There are two major types of methods for detecting AF, first, by calculat‐ ing the R‐R interval (RRI), which represents the heart rate. Majority of the methods in the lit‐ erature are based on RRI [60–68]. Another approach is by using RRI in combination with atrial activity (AA) analysis, which uses heart rate along with atrial ECG pattern. However, in rate controlled patients, the heart rate is externally controlled, in which case, the algorithms that are based on rate will fail to detect AF [69, 70]. Currently available methods working based

have high values in the left atrium [50]. However, there are some studies suggesting that DF

The workstations available today in the EP labs consists of all the components required to assist the clinicians during an ablation procedure, such as the electroanatomic mapping soft‐ ware, mapping and ablation catheters, imaging devices such as fluoroscopy and ultrasound, and other associated software such as the FIRM mapping, DF mapping, and CFAE mapping. Therefore, EP workstation is the best example of a clinical DSS, without which the ablation procedures is nearly impossible to perform today. This demonstrates the significance of clini‐

As discussed earlier, clinical DSSs play an important role in AF ablation treatment. This sec‐ tion presents a detailed discussion of some techniques and algorithms developed by our group that could be potentially implemented in a clinical DSS to assist a clinician for manag‐ ing AF patients. A schematic diagram of the various EP systems along with the integration of

The above three tasks quite covers the necessary tasks required for AF management today. All the algorithms mentioned below were implemented in MATLAB (The MathWorks, Inc.,

Rate control and rhythm control are two strategies for treatment of AF. The goal of rhythm control is to convert an arrhythmic condition to NSR, and that of rate control is to slow down the heart rate when it increases due to AF. Catheter ablation therapy is an example of rhythm control method. Although catheter ablation is a commonly preferred treatment, at an initial stage of AF, doctors recommend rate control therapy. It could be performed using rate control

Detection of AF is a crucial operation for providing timely and appropriate treatment for the patient. AF detection techniques mostly involve application of signal processing algorithms to ECG signals [59]. There are two major types of methods for detecting AF, first, by calculat‐ ing the R‐R interval (RRI), which represents the heart rate. Majority of the methods in the lit‐ erature are based on RRI [60–68]. Another approach is by using RRI in combination with atrial activity (AA) analysis, which uses heart rate along with atrial ECG pattern. However, in rate controlled patients, the heart rate is externally controlled, in which case, the algorithms that are based on rate will fail to detect AF [69, 70]. Currently available methods working based

cal DSS not only in AF, but also in general clinical cardiac electrophysiology.

The role of clinical DSS in AF could be categorized into three important tasks:

**3.1. Developed methods of clinical decision making for AF**

one of the developed clinical DSS system is illustrated in **Figure 8**.

• Identification of successful ablation targets from EGM signals.

*3.1.1. Detection of AF from surface ECG for AF therapy management*

Natick, MA, USA).

drugs or cardiac pacemakers.

• Detection of AF from the surface ECG for AF therapy management.

• Prediction of successful postcardioversion patients with persistent AF.

is not an effective method to analyze AF [37].

104 Computer-aided Technologies - Applications in Engineering and Medicine

**Figure 8.** Schematic of a clinical DSS for AF: The figure shows a schematic of the different components of a typical EP laboratory, integrated with the developed clinical DSS algorithm. The computer screen of the developed DSS illustrates an example of ablation target localization, where the asterisk indicates the estimated rotor location as would be produced by the algorithm described in Section 3.1.3.

on AA exclusively, do not perform well [71]. AA analysis is more complicated in rate‐con‐ trolled patients and could result in false diagnosis mainly because of the low voltage nature of the AA signals compared to the ventricle activities. The idea that the absence of P‐waves is a major indicator of AF is used in most of the AF detection algorithms. However, ECG is a nonstationary signal and P‐waves have relatively less amplitude, so the algorithms using P‐ wave detection encounter many challenges. In this section, we discuss a robust algorithm that we developed for automatic AF detection. The algorithm is independent of heart rate, and is suitable for rate‐controlled patients with paroxysmal AF. The algorithm is capable of fast and reliable detection of AF by using minimum cardiac cycles (heart beats).

The algorithm essentially uses short ECG signals, and applies a feature extraction and machine learning approach to detect whether or not an abnormality exists in the signal. A block diagram of the developed system is shown in **Figure 9**.

The input signal is first subjected to preprocessing. The preprocessing step consists of a third‐ order Butterworth bandpass filter in order to reduce noise. The poles are designed to be at 0.5 and 50 Hz. For each heartbeat, nine features are extracted out of which three are statistical features and the rest are morphological features. The statistical features include variance, skewness, and kurtosis. The morphological features are calculated as the average amplitude

**Figure 9.** AF detection method: Figure shows the block diagram of the developed method of AF detection. PWA = P‐ wave absence, M‐distance = Mahalanobis distance. A, Feature extraction and training; B, classification; C, majority voter post‐processing.

of the P‐wave at different intervals (i.e., first 20 ms or last 10, 20, 20, or 30 ms). A statistical model using a multivariate Gaussian mixture model (GMM) is created using the distribution of the NSR P‐wave features. To generate the training model, the expectation‐maximization (EM) algorithm is iteratively applied [72]. For the testing phase, the same preprocessing is performed and the same features are extracted, but now a majority vote of seven heart beats is employed in order to classify whether the test signal represents An NSR or AF. The algo‐ rithm is evaluated on the MIT‐BIH AF Database from Physiobank [73]. The dataset includes 25 long‐term (10 h) ECG recordings with AF (23 paroxysmal and 2 persistent) and contains 299 AF episodes (about 93.4 h). The proposed algorithm demonstrated that it is capable of classifying AF with a very high sensitivity of 98.08%. The sensitivity was highest among the compared algorithms (see **Table 1**). The specificity was 91.66% which had the second lowest classification error comparatively, which is a significant performance considering that the other methods are RRI‐based approaches.

#### *3.1.2. Prediction of successful postcardioversion patients with persistent AF*

As discussed earlier, DCE cardioversion is one of the treatment options available, in which an electric shock is applied to the chest at certain intervals in order to restore NSR. This treat‐ ment however does not have a high long‐term success rate, although it has more than 90% procedure success. A more well‐defined method to evaluate the effectiveness, and hence

Unreported performance measures are displayed as a dash. Se = Sensitivity, Sp = Specificity, PPV = Positive predictive value, Err = Error rate.

\* Additional algorithm evaluation completed by Larburu et al. is shown on the second row. Optimal window length (WL) was determined by Larburu et al.

† Results include additional processing for the purpose of ectopic beat filtering.

of the P‐wave at different intervals (i.e., first 20 ms or last 10, 20, 20, or 30 ms). A statistical model using a multivariate Gaussian mixture model (GMM) is created using the distribution of the NSR P‐wave features. To generate the training model, the expectation‐maximization (EM) algorithm is iteratively applied [72]. For the testing phase, the same preprocessing is performed and the same features are extracted, but now a majority vote of seven heart beats is employed in order to classify whether the test signal represents An NSR or AF. The algo‐ rithm is evaluated on the MIT‐BIH AF Database from Physiobank [73]. The dataset includes 25 long‐term (10 h) ECG recordings with AF (23 paroxysmal and 2 persistent) and contains 299 AF episodes (about 93.4 h). The proposed algorithm demonstrated that it is capable of classifying AF with a very high sensitivity of 98.08%. The sensitivity was highest among the compared algorithms (see **Table 1**). The specificity was 91.66% which had the second lowest classification error comparatively, which is a significant performance considering that the

**Figure 9.** AF detection method: Figure shows the block diagram of the developed method of AF detection. PWA = P‐ wave absence, M‐distance = Mahalanobis distance. A, Feature extraction and training; B, classification; C, majority voter

As discussed earlier, DCE cardioversion is one of the treatment options available, in which an electric shock is applied to the chest at certain intervals in order to restore NSR. This treat‐ ment however does not have a high long‐term success rate, although it has more than 90% procedure success. A more well‐defined method to evaluate the effectiveness, and hence

other methods are RRI‐based approaches.

106 Computer-aided Technologies - Applications in Engineering and Medicine

post‐processing.

*3.1.2. Prediction of successful postcardioversion patients with persistent AF*

**Table 1.** Comparison of AF detection algorithms including additional evaluation performed by Larburu et al. [74].

predict the likelihood of maintaining NSR after cardioversion for patients who are about to undergo cardioversion, is important, so that the risks involved in cardioversion can be assessed [75]. If such a prediction can be performed, then the doctor could make a decision accordingly on whether to give the patient a cardioversion or not.

This section presents a novel computational method to perform signal analysis of AF patient ECG recorded before cardioversion and predict what the success of the treatment would be. A block diagram of the developed system is shown in **Figure 10**.

The method comprises preprocessing, extraction of features, classification, and then valida‐ tion. The preprocessing is to remove the noise and the isoelectric line perturbations, for which a bandpass filter with cut off frequencies of 0.01 and 50Hz were applied. Following this, AA is extracted from the QRST complex using an average beat subtraction technique. There are a total of seven features extracted from the AA signals. The extracted features are calcu‐ lated from the matching pursuit (MP) time‐frequency representation of the AA signals [76]. A classifier based on quadratic discriminative analysis (QDA) is trained using the extracted features from the training data. The trained classifier was evaluated on an ECG data [77]

**Figure 10.** Prediction of cardioversion outcome: Figure shows the block diagram of the developed prediction method.

with 40 persistent AF patients who had a successful external DCE cardioversion therapy, but after 2‐week follow‐up, 20 patients had maintained NSR (AF‐free) and 20 had a relapse of AF (AF‐relapse). A leave‐one‐out cross‐validation classification produced an outcome of 100% sensitivity and 95% specificity. Compared to the relevant methods in the literature (see **Table 2**), the proposed method resulted in the highest sensitivity and specificity values in predicting the success of DCE cardioversion in AF patients.

#### *3.1.3. Identification of successful ablation targets from EGM signals*

Catheter ablation is the most commonly preferred treatment for AF. However, the success rate of the procedure is only suboptimal due to recurrence of AF especially with persistent type patients. The arrhythmogenic rotors existing outside the PVs could be a potential target for ablation. However, in today's clinical practice there is no well‐defined method to find the location of these rotors. Hence, a system that could estimate the rotor location and hence guide the clinician to the targets would be a significant clinical DSS for AF ablation in today's clinical scenario.

**Table 2.** Comparison of signal processing methods.

The developed method is a novel probabilistic approach to localize rotors [90]. The relevant methods available in the literature for locating rotors involve deterministic, noniterative solutions using time‐domain and frequency‐domain characteristics [54–58]. The developed method uses a Bayesian filtering approach to search for the rotors by iteratively guiding an MPDC toward the center of a rotor. The algorithm was evaluated using a 2D realistic atrial tis‐ sue simulation. We developed a numerical simulation of a 10 cm × 10 cm 2D atrial tissue with a spatial resolution of 0.025 cm and sampling frequency of 500 Hz, using the Nygren human atrial cell model [91]. The numerical simulation was written using FORTRAN programming, and the visualization of the output was implemented in MATLAB. A single stable rotor was initiated and a standard 20 pole Lasso catheter (Biosense Webster, Diamond Bar, CA, USA) was simulated, from which the bipolar EGMs were obtained.

with 40 persistent AF patients who had a successful external DCE cardioversion therapy, but after 2‐week follow‐up, 20 patients had maintained NSR (AF‐free) and 20 had a relapse of AF (AF‐relapse). A leave‐one‐out cross‐validation classification produced an outcome of 100% sensitivity and 95% specificity. Compared to the relevant methods in the literature (see **Table 2**), the proposed method resulted in the highest sensitivity and specificity values in

**Figure 10.** Prediction of cardioversion outcome: Figure shows the block diagram of the developed prediction method.

Catheter ablation is the most commonly preferred treatment for AF. However, the success rate of the procedure is only suboptimal due to recurrence of AF especially with persistent type patients. The arrhythmogenic rotors existing outside the PVs could be a potential target for ablation. However, in today's clinical practice there is no well‐defined method to find the location of these rotors. Hence, a system that could estimate the rotor location and hence guide the clinician to the

targets would be a significant clinical DSS for AF ablation in today's clinical scenario.

**Method Study size Significance Sensitivity Specificity** P‐wave duration, 1997 [78] 35 0.001 73% 71% Heart rate variability, 2001 [79] 93 – 76% 90% Fibrillatory rate, 2003 [80] 44 0.021 – – Clustering of RR intervals, 2004 [81] 66 0.034 – – P‐wave duration, 2005 [82] 118 0.0001 72% 77% P‐wave duration, 2006 [83] 122,493 0.02 90% 21% Fibrillatory rate, 2006 [84] 175 0.0001 79% 80% Fibrillatory rate, 2006 [85] 54 0.002 – – Harmonic decay, 2006 [85] 54 0.0004 92% 47% Sample entropy, 2006 [86] 66 0.02 – – Wavelet transform, 2007 [87] 30 – 100% 89% P‐wave dispersion, 2011 [88] 26 0.001 86% 95% Wavelet sample entropy, 2008 [89] 40 – 95% 93% Proposed MP‐based method, 2014 40 0.005 100% 95%

predicting the success of DCE cardioversion in AF patients.

108 Computer-aided Technologies - Applications in Engineering and Medicine

**Table 2.** Comparison of signal processing methods.

*3.1.3. Identification of successful ablation targets from EGM signals*

As shown in **Figure 11**, the method uses the discrete search space and some EGM characteris‐ tics as input, and gives the probability distribution of the presence of rotor in the search space, from which the rotor location can be estimated as the output. The 2D tissue is the bounded search space and it is discretized to form a grid, the coordinates of which is the first input to the algorithm. The center of the Lasso can be placed anywhere on this grid space. When a catheter is placed at a location, the following EGM characteristics are computed, which is the second input: the first activated bipole (FAB) and the rotor propagation ratio (RPR). FAB denotes the label number of the bipolar electrode of the Lasso that gets activated earliest by the rotor's wavefront. It is mathematically represented by Φ. RPR, represented by Φ, is the ratio of the two following characteristics: total conduction delay and cycle length, the defini‐ tion of which can be found in the publication of our preliminary studies using both simulation and real clinical data [54].

The overall proposed idea is that the clinician places the catheter at any arbitrary endocar‐ dium location and records the EGMs. Then the system calculates FAB and RPR from these EGMs. The search space is discretized and now, both the inputs to the Bayesian algorithm are kept ready. The system first checks if the current catheter location is already a rotor cen‐ ter or not. This is performed using the EGM characteristics and the mathematical criterion, τ≃1, that is established for rotor convergence in our previous work [54]. If it is already a

**Figure 11.** Overview of the guidance algorithm—The inputs are the search space and the two EGM characteristics, RPR (τ) and FAB (Φ) and the output is the estimated location of rotor.

rotor center, then, the clinician is advised to perform ablation, however, if it is not a rotor, i.e., if the convergence condition is not satisfied, then the inputs are provided to the devel‐ oped algorithm to estimate a high probability location. The clinician then moves the catheter to this new location. This completes one iteration of the catheter‐guidance algorithm. The EGMs are recorded again at the new location, and the convergence condition is checked again. If a rotor center is still not found, the new EGMs are processed by the algorithm to output a new estimate, and the catheter is moved by the clinician to this location. This pro‐ cess is repeated until a rotor is located, i.e., until the convergence condition is satisfied. The results of the developed algorithm demonstrate that the time taken to reach convergence is minimized during every iteration. This is mainly because of the updating steps of the recur‐ sive Bayesian filter, which tends to decrease the variance of the posterior distribution with increasing iterations.

The Bayesian algorithm uses the traditional Bayes theorem formulation which states that the posterior is proportional to the product of the prior and the likelihood. The posterior prob‐ ability is nothing but the probability distribution of the presence of rotor in the search space. Here, the task is to determine the likelihood, and for that the EGM characteristics are used. The likelihood in our algorithm is defined as a function of two separate likelihood functions which are each a function of the features calculated from the local EGM recording. The prior is considered to be uniform distribution, meaning that the rotor could be present anywhere on the search space initially. However, in future, the prior could be tuned to any particu‐ lar distribution according to clinicians' expertise and the information from cardiac anatomy. Finally, the algorithm uses the maximum *a posteriori* technique on the posterior distribution to determine the estimated location of a rotor center. This is the procedure during the first iteration; during the second one, the posterior from the previous step becomes the new prior in the Bayesian formulation, and the estimation is carried out again. This process continues until we converge.

Such an iterative framework makes the guidance to be adaptive, which could be much useful to tackle the dynamic nature of rotors. The results of our algorithm are highly promising. The algorithm was tested by placing the catheter in all the possible grid locations on the 2D simu‐ lated tissue. The algorithm achieved 100% convergence, and it took 3.37 ± 1.05 (mean ± sd) number of steps to find a rotor center. In order to assess the usefulness of the developed method, we compared it to a possible method for searching for rotors that could be used in current clinical practice using a Lasso MPDC. It is an unguided approach carried out by ran‐ domly placing the catheter at different locations of the atrium, and checking for convergence at each location. After a couple of placements, if the rotor is not found, the clinician might simply abort the manual search. We implemented this strategy on the simulation, with 20 iterations as the threshold. The convergence in the unguided case was only 34% and the num‐ ber of steps, mean ± sd was 6.58 ± 3.72. This demonstrates the significance of the developed guidance strategy for implementation in a clinical setting.

Some examples of the catheter guidance implemented in MATLAB, is shown in **Figure 12**. The initial arbitrary location of the catheter is s1, and the catheter is advanced (as "dictated" by the algorithm) through s2, and finally converges at s3, which happen to be the rotor center.

**Figure 12.** Some examples of the catheter‐guidance paths using the proposed algorithm on the simulation—Each path is indicated from the first step, s1, to the rotor‐convergence location in s3.

The algorithm's robustness can be clearly seen in the figure, with two major observations: First, the adaptive nature from the distance of the jump; the jump from s1 to s2 is always longer than the subsequent jumps, which means that the method tries to take shorter steps when it "feels" that it is getting close to a rotor center. The second observation is the adaptiv‐ ity with respect to the direction; the direction of the path diverts toward the center of the rotor with increasing iterations. Additionally, it is also seen that the convergence points (red dots in the figure) are very close to the center (according to the developed simulation). The devel‐ oped algorithm however is only a preliminary design and has to be improvised in future to tackle practical scenarios such as heterogeneous wave propagation, rotors with different cycle lengths, rotors with meandering center, etc.

## **4. Conclusion**

rotor center, then, the clinician is advised to perform ablation, however, if it is not a rotor, i.e., if the convergence condition is not satisfied, then the inputs are provided to the devel‐ oped algorithm to estimate a high probability location. The clinician then moves the catheter to this new location. This completes one iteration of the catheter‐guidance algorithm. The EGMs are recorded again at the new location, and the convergence condition is checked again. If a rotor center is still not found, the new EGMs are processed by the algorithm to output a new estimate, and the catheter is moved by the clinician to this location. This pro‐ cess is repeated until a rotor is located, i.e., until the convergence condition is satisfied. The results of the developed algorithm demonstrate that the time taken to reach convergence is minimized during every iteration. This is mainly because of the updating steps of the recur‐ sive Bayesian filter, which tends to decrease the variance of the posterior distribution with

110 Computer-aided Technologies - Applications in Engineering and Medicine

The Bayesian algorithm uses the traditional Bayes theorem formulation which states that the posterior is proportional to the product of the prior and the likelihood. The posterior prob‐ ability is nothing but the probability distribution of the presence of rotor in the search space. Here, the task is to determine the likelihood, and for that the EGM characteristics are used. The likelihood in our algorithm is defined as a function of two separate likelihood functions which are each a function of the features calculated from the local EGM recording. The prior is considered to be uniform distribution, meaning that the rotor could be present anywhere on the search space initially. However, in future, the prior could be tuned to any particu‐ lar distribution according to clinicians' expertise and the information from cardiac anatomy. Finally, the algorithm uses the maximum *a posteriori* technique on the posterior distribution to determine the estimated location of a rotor center. This is the procedure during the first iteration; during the second one, the posterior from the previous step becomes the new prior in the Bayesian formulation, and the estimation is carried out again. This process continues

Such an iterative framework makes the guidance to be adaptive, which could be much useful to tackle the dynamic nature of rotors. The results of our algorithm are highly promising. The algorithm was tested by placing the catheter in all the possible grid locations on the 2D simu‐ lated tissue. The algorithm achieved 100% convergence, and it took 3.37 ± 1.05 (mean ± sd) number of steps to find a rotor center. In order to assess the usefulness of the developed method, we compared it to a possible method for searching for rotors that could be used in current clinical practice using a Lasso MPDC. It is an unguided approach carried out by ran‐ domly placing the catheter at different locations of the atrium, and checking for convergence at each location. After a couple of placements, if the rotor is not found, the clinician might simply abort the manual search. We implemented this strategy on the simulation, with 20 iterations as the threshold. The convergence in the unguided case was only 34% and the num‐ ber of steps, mean ± sd was 6.58 ± 3.72. This demonstrates the significance of the developed

Some examples of the catheter guidance implemented in MATLAB, is shown in **Figure 12**. The initial arbitrary location of the catheter is s1, and the catheter is advanced (as "dictated" by the algorithm) through s2, and finally converges at s3, which happen to be the rotor center.

guidance strategy for implementation in a clinical setting.

increasing iterations.

until we converge.

This chapter reviewed some of the clinical decision support systems that have been developed and used in the therapy management of AF patient. Specifically, we discussed methods for detection of AF episodes from surface ECG signals, prediction of the success of DCE cardio‐ version from ECG signals, and identification of ablation targets for AF catheter ablation. The clinical DSS development is taking place in the right pace, the research field is witnessing promising methods that could turn into reliable clinical DSS technologies to assist the clini‐ cians in the process of AF treatment and management in order to enhance the health of mil‐ lions of patients afflicted by this debilitating rhythm disorder so that they can live longer and more fulfilling lives.

## **Acknowledgements**

A part of the work was funded by the National Institutes of Health R15HL127663 grant to Dr. Ghoraani.

## **Conflict of interest**

The authors declare no conflict of interests.

## **Author details**

Prasanth Ganesan1 , Mark Sterling<sup>2</sup> , Steven Ladavich3 and Behnaz Ghoraani<sup>1</sup> \*

\*Address all correspondence to<EMAIL_ADDRESS>
1 Department of Computer and Electrical Engineering, Florida Atlantic University, Boca Raton, Florida, USA

2 Department of Computer Science, Nazarbayev University, Astana, Republic of Kazakhstan

3 Department of Biomedical Engineering, Rochester Institute of Technology, Rochester, New York, USA

## **References**

[7] Power DJ. Decision support systems: a historical overview. In: Handbook on Decision Support Systems 1. Springer Berlin Heidelberg; 2008. pp. 121–140.

**Acknowledgements**

**Conflict of interest**

**Author details**

Prasanth Ganesan1

Raton, Florida, USA

New York, USA

**References**

The authors declare no conflict of interests.

, Mark Sterling<sup>2</sup>

112 Computer-aided Technologies - Applications in Engineering and Medicine

\*Address all correspondence to<EMAIL_ADDRESS>
Hall PTR Upper Saddle River, NJ, USA; 1990.

Management Science. 1966;12(8):B‐363–B‐381.

Management Science. 1967;13(6):B‐342–B‐358.

Norton & Company; 2001. pp. 64–90.

Dr. Ghoraani.

A part of the work was funded by the National Institutes of Health R15HL127663 grant to

, Steven Ladavich3

1 Department of Computer and Electrical Engineering, Florida Atlantic University, Boca

2 Department of Computer Science, Nazarbayev University, Astana, Republic of Kazakhstan 3 Department of Biomedical Engineering, Rochester Institute of Technology, Rochester,

[1] Turban E. Decision Support and Expert Systems: Management Support Systems. Prentice

[2] Raymond RC. Use of the time‐sharing computer in business planning and budgeting.

[3] Turban E. The use of mathematical models in plant maintenance decision making.

[4] Shortliffe E. Computer‐Based Medical Consultations: MYCIN. vol. 2. American Elsevier

[5] Ancelin C, Le P, De Saint‐Quentin S, Villatte N. EXPRESS: an expert system to per‐ form system safety studies. In: Artificial Intelligence and Other Innovative Computer

[6] Engelbart DC. Augmenting human intellect: a conceptual framework (1962). PACKER, Randall and JORDAN, Ken Multimedia from Wagner to Virtual Reality. New York: WW

Publishing Company, Inc. 52 Vanderbilt Avenue, New York, NY 10017; 2012.

Applications in the Nuclear Industry. Springer US; 1988. pp. 761–768.

and Behnaz Ghoraani<sup>1</sup>

\*

Ablation of Atrial Fibrillation. Developed in partnership with the European Heart Rhythm Association (EHRA), a registered branch of the European Society of Cardiology (ESC) and the European Cardiac Arrhythmia Society (ECAS); and in collaboration with the American College of Cardiology (ACC), American Heart Association (AHA), the Asia Pacific Heart Rhythm Society (APHRS), and the Society of Thoracic Surgeons (STS). Endorsed by the governing bodies of the American College of Cardiology Foundation, the American Heart Association, the European Cardiac Arrhythmia Society, the European Heart Rhythm Association, the Society of Thoracic Surgeons, the Asia Pacific Heart Rhythm Society, and the Heart Rhythm Society. Europace. 2012;14(4):528–606.

[23] Doi K. Current status and future potential of computer‐aided diagnosis in medical imag‐

[24] Doi K. Computer‐aided diagnosis in medical imaging: historical review, current status and future potential. Computerized Medical Imaging and Graphics. 2007;31(4):198–211.

[25] Acharya UR, Sudarshan VK, Adeli H, Santhosh J, Koh JEW, Adeli A. Computer‐aided diagnosis of depression using EEG signals. European Neurology. 2015;73(5–6):329–336.

[26] Huynh B, Drukker K, Giger M. MO‐DE‐207B‐06: Computer‐aided diagnosis of breast ultrasound images using transfer learning from deep convolutional neural networks.

[27] Hirschauer TJ, Adeli H, Buford JA. Computer‐aided diagnosis of Parkinsons disease using enhanced probabilistic neural network. Journal of Medical Systems. 2015;39(11):1–12.

[28] Mozaffarian D, Benjamin EJ, Go AS, Arnett DK, Blaha MJ, Cushman M, et al. Heart disease and stroke statistics update a report from the American Heart Association.

[29] Chugh SS, Roth GA, Gillum RF, Mensah GA. Global burden of atrial fibrillation in devel‐

[30] Issa Z, Miller JM, Zipes DP. Clinical Arrhythmology and Electrophysiology: A Companion to Braunwald's Heart Disease: Expert Consult: online and print. Elsevier Saunders, 1600 John F. Kennedy Bld. Ste 1800, Philadelphia, PA 19103‐2899; 2012. [31] Hurst JW. Naming of the waves in the ECG, with a brief account of their genesis.

[32] Hlaing T, DiMino T, Kowey PR, Yan G. ECG repolarization waves: their genesis and clinical implications. Annals of Noninvasive Electrocardiology. 2005;10(2):211–223. [33] Ganesan P. Characterization of cardiac electrogram signals during atrial fibrillation. Thesis, Rochester Institute of Technology. 2015; Accessed fromhttp://scholarworks.rit.

[34] Atkielski A, Wikimedia. Sinus Rhythm Labels; 2007. Available from: https://upload.wiki‐

[35] Houghton A, Gray D. Making Sense of the ECG: A Hands‐On Guide. CRC Press, Taylor & Francis Group, 6000 Broken Sound Parkway NW, Duite 300, Boca Raton, FL 33487–

[37] Calkins H, Kuck KH, Cappato R, Brugada J, Camm AJ, Chen SA, et al. 2012 HRS/ EHRA/ECAS Expert Consensus Statement on Catheter and Surgical Ablation of Atrial Fibrillation: Recommendations for Patient Selection, Procedural Techniques, Patient Management and Follow‐up, Definitions, Endpoints, and Research Trial Design: A report of the Heart Rhythm Society (HRS) Task Force on Catheter and Surgical

[36] Fenton FH, Cherry EM, Glass L. Cardiac arrhythmia. Scholarpedia. 2008;3(7):1665.

oped and developing nations. Global Heart. 2014;9(1):113–119.

media.org/wikipedia/commons/9/9e/SinusRhythmLabels.svg

ing. The British Journal of Radiology. 2014; 78, pages s3–s19.

Medical Physics. 2016;43(6):3705–3705.

114 Computer-aided Technologies - Applications in Engineering and Medicine

Circulation. 2015;131(4):e29–e322.

Circulation. 1998;98(18):1937–1942.

edu/theses/8760

2742; 2014.

[61] Cerutti S, Mainardi LT, Porta A, Bianchi AM. Analysis of the dynamics of RR interval series for the detection of atrial fibrillation episodes. In: Computers in Cardiology. 1997;24, pp. 77–80.

[49] Narayan SM, Krummen DE, Shivkumar K, Clopton P, Rappel WJ, Miller JM. Treatment of atrial fibrillation by the ablation of localized sources: CONFIRM (conventional abla‐ tion for atrial fibrillation with or without focal impulse and rotor modulation) trial.

[50] Pandit SV, Jalife J. Rotors and the dynamics of cardiac fibrillation. Circulation Research.

[51] Nash MP, Mourad A, Clayton RH, Sutton PM, Bradley CP, Hayward M, et al. Evidence for multiple mechanisms in human ventricular fibrillation. Circulation.

[52] Bray MA, Wikswo JP. Considerations in phase plane analysis for nonstationary reen‐

[53] Baykaner T, Lalani GG, Schricker A, Krummen DE, Narayan SM. Mapping and ablat‐ ing stable sources for atrial fibrillation: summary of the literature on focal impulse and rotor modulation (FIRM). Journal of Interventional Cardiac Electrophysiology.

[54] Ganesan P, Cherry EM, Pertsov AM, Ghoraani B. Characterization of electrograms from multipolar diagnostic catheters during atrial fibrillation. BioMed Research International.

[55] Arunachalam SP, Annoni EM, Mulpuru SK, Friedman PA, Tolkacheva EG. Novel mul‐ tiscale frequency approach to identify the pivot point of the rotor. Journal of Medical

[56] Shariat MH, Gazor S, Redfearn D. Computationally efficient method for localizing the spiral rotor source using synthetic intracardiac electrograms during atrial fibrillation. In: Engineering in Medicine and Biology Society (EMBC), 2015 37th Annual International

[57] Orozco‐Duque A, Duque SI, Ugarte JP, Tobon C, Novak D, Kremen V, et al. Fractionated electrograms and rotors detection in chronic atrial fibrillation using model‐based cluster‐ ing. In: 2014 36th Annual International Conference of the IEEE Engineering in Medicine

[58] Balasundaram K, Umapathy K, Jeyaratnam J, Niri A, Masse S, Farid T, et al. Tracking rotors with minimal electrodes: Modulation Index Based Strategy. Circulation:

[59] Langley P, Dewhurst M, Di Marco L, Adams P, Dewhurst F, Mwita J, et al. Accuracy of algorithms for detection of atrial fibrillation from short duration beat interval record‐

[60] Moody GB, Mark RG. A new method for detecting atrial fibrillation using RR intervals.

Journal of the American College of Cardiology. 2012;60(7):628–636.

116 Computer-aided Technologies - Applications in Engineering and Medicine

trant cardiac behavior. Physical Review E. 2002;65(5):051902.

2013;112(5):849–862.

2006;114(6):536–542.

2014;40:237–244.

Elsevier, US 2015;2015.

Devices. 2016;10(2):020948.

Conference of the IEEE. IEEE; 2015. pp. 4483–4486.

and Biology Society. IEEE; 2014. pp. 1579–1582.

Computers in Cardiology. 1983;10(1):227–230.

Arrhythmia and Electrophysiology; 2015; 8, pp. 447–455.

ings. Medical Engineering and Physics. 2012;34(10):1441–1447.

[87] Watson JN, Addison PS, Uchaipichat N, Shah AS, Grubb NR. Wavelet transform analy‐ sis predicts outcome of DC cardioversion for atrial fibrillation patients. Computers in Biology and Medicine. 2007;37(4):517–523.

[75] Sterling M, Huang DT, Ghoraani B. Developing a new computer‐aided clinical decision support system for prediction of successful postcardioversion patients with persistent atrial fibrillation. Computational and Mathematical Methods in Medicine. vol. 2015,

[76] Mallat SG, Zhang Z. Matching pursuits with time‐frequency dictionaries. IEEE

[77] Couderc JP. The telemetric and Holter ECG warehouse initiative (THEW): a data reposi‐ tory for the design, implementation and validation of ECG‐related technologies. In: Conference proceedings: Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE Engineering in Medicine and Biology Society.

[78] Opolski G, Scislo P, Stanislawska J, Górecki A, Steckiewicz R, Torbicki A. Detection of patients at risk for recurrence of atrial fibrillation after successful electrical cardioversion by signal‐averaged P‐wave ECG. International Journal of Cardiology. 1997;60(2):181–185.

[79] Lombardi F, Colombo A, Basilico B, Ravaglia R, Garbin M, Vergani D, et al. Heart rate variability and early recurrence of atrial fibrillation after electrical cardioversion. Journal

[80] Bollmann A, Husser D, Steinert R, Stridh M, Soernmo L, Olsson SB, et al. Echocardiographic and electrocardiographic predictors for atrial fibrillation recurrence following cardio‐

[81] Van Den Berg MP, Van Noord T, Brouwer J, Haaksma J, Van Veldhuisen DJ, Crijns HJ, et al. Clustering of RR intervals predicts effective electrical cardioversion for atrial fibril‐

[82] Budeus M, Hennersdorf M, Perings C, Wieneke H, Erbel R, Sack S. Prediction of the recurrence of atrial fibrillation after successful cardioversion with P wave signal‐aver‐

[83] Raitt MH, Volgman AS, Zoble RG, Charbonneau L, Padder FA, O'Hara GE, et al. Prediction of the recurrence of atrial fibrillation after cardioversion in the atrial fibrilla‐ tion follow‐up investigation of rhythm management (AFFIRM) study. American Heart

[84] Holmqvist F, Stridh M, Waktare JE, Sörnmo L, Olsson SB, Meurling CJ. Atrial fibrillatory rate and sinus rhythm maintenance in patients undergoing cardioversion of persistent

[85] Holmqvist F, Stridh M, Waktare JE, Roijer A, Sörnmo L, Platonov PG, et al. Atrial fibril‐ lation signal organization predicts sinus rhythm maintenance in patients undergoing

[86] Petersson R, Sandberg F, Platonov PG, Holmqvist F. Noninvasive estimation of orga‐ nization in atrial fibrillation as a predictor of sinus rhythm maintenance. Journal of

version. Journal of Cardiovascular Electrophysiology. 2003;14(s10):S162–S165.

lation. Journal of Cardiovascular Electrophysiology. 2004;15(9):1027–1033.

aged ECG. Annals of Noninvasive Electrocardiology. 2005;10(4):414–419.

atrial fibrillation. European Heart Journal. 2006;27(18):2201–2207.

cardioversion of atrial fibrillation. Europace. 2006;8(8):559–565.

Article ID 527815, 10 pages, doi:10.1155/2015/527815

118 Computer-aided Technologies - Applications in Engineering and Medicine

Transactions on Signal Processing. 1993;41(12):3397–3415.

Conference. vol. 2010. NIH Public Access; 2010. p. 6252.

of the American College of Cardiology. 2001;37(1):157–162.

Journal. 2006;151(2):390–396.

Electrocardiology. 2011;44(2):171–175.

**Chapter 6 Provisional chapter**

#### **Computer-Aided Biosensor Design Computer-Aided Biosensor Design**

Yu-Chen Lo, Ren Gui, Hiroshi Honda and Jorge Z. Torres Yu-Chen Lo, Ren Gui, Hiroshi Honda and Jorge Z. Torres

Additional information is available at the end of the chapter Additional information is available at the end of the chapter

http://dx.doi.org/5772/65038

#### **Abstract**

Amperometric biosensors are widely used in point-of-care medical devices that help patients control blood glucose and cholesterol levels in an effective and convenient way. On the other hand, computer-aided technologies for biosensor design remain an actively developing field. In this chapter, we present a computational model for biosensor design that uses a reaction-diffusion equation. We have successfully applied this model to simulate cholesterol analysis based on a multienzyme system. Furthermore, we show that this computer-aided approach can be used to optimize biosensor performance. This model can be applied to industry-grade biosensor development and can be easily extended to model multiple types of biosensors for a wide array of clinical applications.

**Keywords:** biosensor, multienzyme system, reaction-diffusion equation, computeraided modeling and simulation

## **1. Introduction**

High blood cholesterol is a major risk factor of cardiovascular disease including hypertension, hyperthyroidism, anemia, and coronary artery disease. Many prevention strategies including lifestyle modification, diet, and cholesterol-lowering drugs have been indicated for early prevention and management of the diseases. For adults, the normal cholesterol level is less than 200 mg/dL, and greater than 200 mg/dL is associated with increased risk of disease [1]. For detecting whole blood cholesterol levels and diagnosing cardiovascular disease, a biosensor is a convenient way to monitor cholesterol levels because the device is compact and easy to carry, and the diagnostic results can be obtained in seconds. A biosensor is defined as an analytical device that incorporates a transducer and biologically active material within

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

it [2, 3]. A biosensor obtains a signal that is proportional to the concentration of the analyte measurement. There are many types of biosensors including potentiometric, optical, and amperometric types [4–6]. The most common type is the amperometric biosensor, which measures the concentration of the analytes based on current. To use a biosensor for cholesterol-level determination, the patient uses a lancet to prick a small nick on the fingertip and apply a small drop of blood (around 10–25 μl) on the reactive area of the test strip. The test strip is coated with immobilized enzymes, which initiate the enzymatic reaction necessary for signal production.

## **2. Biosensor design**

Although enzyme composition is optimized for bulk reagent laboratory testing, the composition does not apply to biosensors used for point-of-care applications. As a result, the enzyme composition needs to be reformulated and optimized for a specific biosensor platform. The enzyme formulation step is often labor intensive and time-consuming and requires tremendous amount of manual testing and calibration to achieve optimal performance. The construction of a mathematical model that can model biosensor kinetics and extract important system parameters to improve the biosensor performance will facilitate a more cost-effective biosensor design. The three most important criteria in biosensor characterization are sensitivity, measurement time, and measurement range [7]. Sensitivity measures the ability of the biosensor to differentiate close analytical values with minimal error and can be quantified based on the slope of the current density versus substrate calibration curve. The measurement time measures how fast the biosensor signal reaches steady state and is able to establish a linear relationship between substrate concentration and current magnitude. The measurement range is the range of the substrate concentration in the linear region of the calibration curve and will often include both the normal and pathological values.

To design a biosensor, a mathematical model can be constructed based on the design specification including input (type of analyte and its concentration), signal output (current density), and system parameters (enzyme compositions). The model can be then validated by comparing its signal output with the real-time biosensor data to test its accuracy. Meanwhile, several assumptions of the model including boundary conditions and initial values can be adjusted based on the specification. Once the model is validated, system parameters such as enzyme compositions and reaction area can be further adjusted to optimize biosensor signal.

#### **2.1. Modeling enzymatic reactions**

The chemical reaction of the amperometric biosensor is known as a heterogeneous reaction that occurs at the interface between the membrane and sensing electrode through the immobilized enzymes (**Figure 1**) [8].

As a result, the transient response of the amperometric biosensor can be modeled using a reaction-diffusion equation with the reaction term corresponding to the Michaelis-Menten kinetics of the enzyme reactions [9]:

Computer-Aided Biosensor Design http://dx.doi.org/5772/65038 123

$$\frac{\partial \left[C\right]}{\partial t} = D \mathbf{e} \frac{\partial^2 \left[C\right]}{\partial \mathbf{x}^2} + R\_\mathbf{e} \tag{1}$$

where

it [2, 3]. A biosensor obtains a signal that is proportional to the concentration of the analyte measurement. There are many types of biosensors including potentiometric, optical, and amperometric types [4–6]. The most common type is the amperometric biosensor, which measures the concentration of the analytes based on current. To use a biosensor for cholesterol-level determination, the patient uses a lancet to prick a small nick on the fingertip and apply a small drop of blood (around 10–25 μl) on the reactive area of the test strip. The test strip is coated with immobilized enzymes, which initiate the enzymatic reaction necessary

Although enzyme composition is optimized for bulk reagent laboratory testing, the composition does not apply to biosensors used for point-of-care applications. As a result, the enzyme composition needs to be reformulated and optimized for a specific biosensor platform. The enzyme formulation step is often labor intensive and time-consuming and requires tremendous amount of manual testing and calibration to achieve optimal performance. The construction of a mathematical model that can model biosensor kinetics and extract important system parameters to improve the biosensor performance will facilitate a more cost-effective biosensor design. The three most important criteria in biosensor characterization are sensitivity, measurement time, and measurement range [7]. Sensitivity measures the ability of the biosensor to differentiate close analytical values with minimal error and can be quantified based on the slope of the current density versus substrate calibration curve. The measurement time measures how fast the biosensor signal reaches steady state and is able to establish a linear relationship between substrate concentration and current magnitude. The measurement range is the range of the substrate concentration in the linear region of the calibration curve and will

To design a biosensor, a mathematical model can be constructed based on the design specification including input (type of analyte and its concentration), signal output (current density), and system parameters (enzyme compositions). The model can be then validated by comparing its signal output with the real-time biosensor data to test its accuracy. Meanwhile, several assumptions of the model including boundary conditions and initial values can be adjusted based on the specification. Once the model is validated, system parameters such as enzyme

The chemical reaction of the amperometric biosensor is known as a heterogeneous reaction that occurs at the interface between the membrane and sensing electrode through the

As a result, the transient response of the amperometric biosensor can be modeled using a reaction-diffusion equation with the reaction term corresponding to the Michaelis-Menten

compositions and reaction area can be further adjusted to optimize biosensor signal.

for signal production.

**2. Biosensor design**

often include both the normal and pathological values.

122 Computer-aided Technologies - Applications in Engineering and Medicine

**2.1. Modeling enzymatic reactions**

immobilized enzymes (**Figure 1**) [8].

kinetics of the enzyme reactions [9]:

$$R\_v = \pm \frac{V\_{\text{max}} \left\lfloor C \right\rfloor}{k s + \left\lceil C \right\rceil} \tag{2}$$

Here, *V*max is the maximal rate of the enzymatic reaction, which is proportional to the total amount of enzyme loaded on the membrane.

$$V\_{\text{max}} = k\_{\text{car}} \text{[}Et\text{]}\tag{3}$$

The current density of the biosensor can also be calculated based on the flux of the electroactive species at the sensing membrane:

$$i(t) = zFDp \left. \frac{\partial [P]}{\partial \alpha} \right|\_{\alpha=0} \tag{4}$$

*z* is the number of electrons transferred per species, *D* is the diffusion constant for the electroactive species, and F is the Faraday constant (F = 96485 C/mole).

**Figure 1.** Layout of a biosensor test strip. The test strip consists of three layers: diffusion layer, membrane, and sensing electrode. The serum analytes go from the top layer through the diffusion layer and reach the membrane where enzymatic reactions take place. The sensing electrode converts electrons from the electron-donating species in the reaction to a current signal that can be used to quantify the amount of serum analytes.

#### **2.2. One-enzyme reactions**

The physical significance of a diffusion layer accounts for the serum diffusion. Different blood volumes can produce significant variance in the biosensor signal. As a result, volume control becomes an important issue in biosensor design. The reactions that occur at the membrane region are

$$
\Delta S + E \Leftrightarrow ES \to E + P \tag{5}
$$

$$P \to P' \pm e - \tag{6}$$

In the first reaction, *S* represents the applied analyte, and *E* represents the catalytic enzyme. The enzyme first converts the analyte into an electro-active species *P*. The second reaction is the electron-donating step in which the electro-active species donate electrons at the membrane electrode. Since the electro-reaction is very fast, the electro-reaction term is ignored in the model. The reactions can be captured using the system reaction-diffusion equation in nondimensionalized form:

$$\frac{\partial \left[S\right]}{\partial T} = \frac{\partial^2 \left[S\right]}{\partial X^2} - \phi^2 \frac{\left[S\right]}{1 + \left[S\right]}\tag{7}$$

$$\frac{\partial \left[P\right]}{\partial T} = D \frac{\partial^2 \left[P\right]}{\partial X^2} + \phi^2 \frac{\left[S\right]}{1 + \left[S\right]}\tag{8}$$

where <sup>=</sup> lm , <sup>=</sup> lm , and lm is the membrane thickness. The behavior of the system equation is dictated by the parameter called the Thiele modulus (*ϕ*):

$$\phi = \text{lm}\sqrt{\frac{V\_{\text{max}}}{Dk s}}\tag{9}$$

The Thiele modulus represents the ratio of the diffusion time (1/ ) over reaction time (1/ and is often the reaction characteristic of a particular enzyme in the reaction [10]. For an amperometric biosensor, the Thiele modulus is also directly proportional to the square roots of the enzyme concentration as well as the length of the membrane layer (lm), which is important for reaction control. Another important parameter to characterize an enzymatic diffusion equation is the Biot number:

$$B\_i = \frac{\text{Im} / \text{Dm}}{\text{le} / \text{De}} \tag{10}$$

The Biot number measures the ratio of internal diffusion and external diffusion and can be approximated by the ratio of the thickness between the inner membrane (lm) and outer diffusion layer (le) for a constant diffusion constant D. Previous studies indicate that the biosensor achieves the best sensitivity at high Thiele modulus and high Biot number, that is, under the conditions where the reaction is in internal diffusion control with minimal external diffusion. To minimize external serum diffusion, the pipette blood volume can be controlled based on the geometric design of the reaction region.

## **3. Application**

**2.2. One-enzyme reactions**

124 Computer-aided Technologies - Applications in Engineering and Medicine

dimensionalized form:

where <sup>=</sup>

lm , <sup>=</sup>

diffusion equation is the Biot number:

region are

The physical significance of a diffusion layer accounts for the serum diffusion. Different blood volumes can produce significant variance in the biosensor signal. As a result, volume control becomes an important issue in biosensor design. The reactions that occur at the membrane

'

In the first reaction, *S* represents the applied analyte, and *E* represents the catalytic enzyme. The enzyme first converts the analyte into an electro-active species *P*. The second reaction is the electron-donating step in which the electro-active species donate electrons at the membrane electrode. Since the electro-reaction is very fast, the electro-reaction term is ignored in the model. The reactions can be captured using the system reaction-diffusion equation in non-

[ ] [ ] [ ]

[ ] [ ] [ ]

*TX S* f

> max lm *<sup>V</sup> Dks*

The Thiele modulus represents the ratio of the diffusion time (1/ ) over reaction time (1/ and is often the reaction characteristic of a particular enzyme in the reaction [10]. For an amperometric biosensor, the Thiele modulus is also directly proportional to the square roots of the enzyme concentration as well as the length of the membrane layer (lm), which is important for reaction control. Another important parameter to characterize an enzymatic

*P PS*

2

= +

*D*

f

¶ ¶

is dictated by the parameter called the Thiele modulus (*ϕ*):

2 <sup>2</sup> 1 *SS S TX S* f

> 2 <sup>2</sup> 1

2

[ ]

[ ]

lm , and lm is the membrane thickness. The behavior of the system equation

¶ ¶ = - ¶¶ + (7)

¶¶ + (8)

= (9)

*S E ES E P* +Û ®+ (5)

*P Pe* ® ±- (6)

#### **3.1. Enzymatic determination of serum cholesterol**

Methods for determining serum cholesterol levels have been developed for prepared aqueous reagents (**Figure 2**) [11]. In the first step, the cholesterol ester is converted to free cholesterol and fatty acids using the cholesterol esterase (CE) enzyme. The free cholesterol is then converted to electro-active species hydrogen peroxide and cholest-4-en-3-one using the cholesterol oxidase (CO) enzyme. Hydroxide is an electron-donating species that can be used to generate an electrical current to quantify the amount of cholesterol in the blood serum. The determination scheme is shown below.

**Figure 2.** Enzymatic reaction for the quantification of cholesterol. In the first step, cholesterol ester is converted to free cholesterol and fatty acid through hydrolysis by the cholesterol esterase (CE) enzyme. In the second step, free cholesterol is converted to cholest-4-en-3-one and the electron-donating species hydrogen peroxide by the cholesterol oxidase (CO) enzyme.

#### **3.2. Modeling multienzyme systems**

To model the biosensor kinetics for cholesterol determination, a multienzyme model is constructed based on the one-enzyme model. Assume the biosensor is operated under internal diffusion control, that is, Biot number is infinite, and the serum diffusion layer thickness is small compared with the membrane thickness. The biosensor kinetics is simulated only in the membrane region.

#### *3.2.1. Multienzyme reaction-diffusion equation*

The conversion of cholesterol ester to hydroxide is a two-step enzymatic reaction and thus can be captured using a multienzyme reaction-diffusion equation. In non-dimensionalized form

$$\frac{\partial \left[S1\right]}{\partial T} = \frac{\partial^2 \left[S1\right]}{\partial X^2} - \left(\phi 1\right)^2 \frac{\left[S2\right]}{1 + \left[S2\right]}\tag{11}$$

$$\frac{\partial \left[S2\right]}{\partial T} = \frac{\partial^2 \left[S2\right]}{\partial X^2} - \left(\phi 2\right)^2 \frac{\left[S2\right]}{1 + \left[S2\right]} + \left(\phi 1\right)^2 \frac{\left[S1\right]\left(\frac{Ks}{Ks2}\right)}{1 + \left[S1\right]}\tag{12}$$

$$\frac{\partial \left[P\right]}{\partial T} = D \frac{\partial^2 \left[P\right]}{\partial X^2} + \left(\phi \mathbf{2}\right)^2 \frac{\left[S\mathbf{2}\right]}{1 + \left[S\mathbf{2}\right]}\tag{13}$$

#### *3.2.2. Boundary and initial conditions*

The multienzyme reaction-diffusion equation for cholesterol analysis is a system partial differential equation, and to solve the equation, six boundary conditions and six initial conditions are required.

#### *<IP_ADDRESS>. Boundary conditions 1 and 2*

Since the cholesterol ester and free cholesterol are nonreactive at the electrode surface, their flux is zero:

$$
\left.\frac{\partial\left[S1\right]}{\partial\mathbf{x}}\right|\_{\mathbf{x}=0} = \left.\frac{\partial\left[S2\right]}{\partial\mathbf{x}}\right|\_{\mathbf{x}=0} = \mathbf{0} \tag{14}
$$

#### *<IP_ADDRESS>. Boundary condition 3*

small compared with the membrane thickness. The biosensor kinetics is simulated only in the

The conversion of cholesterol ester to hydroxide is a two-step enzymatic reaction and thus can be captured using a multienzyme reaction-diffusion equation. In non-dimensionalized form

2

[ ] ( )

2

2

The multienzyme reaction-diffusion equation for cholesterol analysis is a system partial differential equation, and to solve the equation, six boundary conditions and six initial

Since the cholesterol ester and free cholesterol are nonreactive at the electrode surface, their

0 0

= = ¶ ¶ = = ¶ ¶ *x x S S*

0

*x x* (14)

[ ] [ ]

1 2

1 2 11

[ ]

2

1 2

 f

2 2

[ ]

[ ]

¶ ¶ = - ¶¶ + (11)

[ ]

¶¶ + (13)

æ ö

(12)

1 2

[ ] [ ] ( ) [ ]

<sup>1</sup> <sup>1</sup> 22 2 <sup>2</sup> 2 1

*Ks <sup>S</sup> SS S Ks*

ç ÷ ¶ ¶ è ø =- + ¶¶ + +

*TX S S*

[ ] [ ] ( ) [ ]

*TX S* f

*PP S*

2

*SS S TX S* f

2 11 2 1

2

[ ] [ ] ( ) [ ]

¶ ¶

f

2

= +

*D*

[S1] = cholesterol ester, [S2] = cholesterol, and [P] = hydrogen peroxide

2

2

ks1 (cholesterol esterase catalytic constant) = 1.9 × 10−5 [M] ks2 (cholesterol oxidase catalytic constant) = 3.5 × 10−4 [M]

membrane region.

*3.2.1. Multienzyme reaction-diffusion equation*

126 Computer-aided Technologies - Applications in Engineering and Medicine

*ϕ*1 = cholesterol esterase Thiele modulus *ϕ*2 = cholesterol oxidase Thiele modulus

*V* max = 30 u\*mg for both enzymes

*3.2.2. Boundary and initial conditions*

*<IP_ADDRESS>. Boundary conditions 1 and 2*

conditions are required.

flux is zero:

Since hydrogen peroxide is consumed very fast, thus there is zero concentration at the electrode surface:

$$\left.P\right|\_{\iota=0} = 0\tag{15}$$

#### *<IP_ADDRESS>. Boundary conditions 4–6*

Assume the concentration of cholesterol ester at the top of the membrane maintains constant concentration and the concentration of cholesterol and hydrogen peroxide are zero:

$$\left.S\_1\right|\_{\mathbf{x}\leftarrow\mathbf{lm}} = S\_O \tag{16}$$

$$\left.S\_2\right|\_{x=\text{lm}} = \mathbf{0} \tag{17}$$

$$\left.P\right|\_{x=\text{lm}} = 0\tag{18}$$

#### *3.2.3. Initial condition*

At the initial time point, there is only cholesterol ester S1 on top of the sensing membrane with a concentration S0. As the reaction has yet to start, there is no cholesterol S2 or hydrogen peroxide P presented in the membrane:

$$\text{At } t = 0 \quad \text{S1} = \text{S0} \quad X \ge 1 \tag{19}$$

$$\text{L'S2} = 0 \quad 0 \le X < 1 \tag{20}$$

$$\mathbf{P} = \mathbf{0} \qquad \forall X \tag{21}$$

#### **3.3. MATLAB simulation**

The biosensor performance for multienzyme reaction cholesterol determination is simulated by solving a system of partial differential equations using MATLAB. The system of partial differential equations is casted into matrix form and solved using the MATLAB function PDEPE by assuming slab geometry. All units are non-dimensionalized. Note that *X* = 1 is the top of the membrane and *X* = 0 is the side of sensing electrode. For the simplified model, the Thiele modulus for both CE and CO is 1.

The solution of the system reaction-diffusion equation for cholesterol analysis is shown in **Figures 3**–**5**. The concentration profile shows that only 25% of the cholesterol ester has been converted to cholesterol at steady-state time and maintains at a constant concentration. On the other hand, the free cholesterol achieves 100% conversion to hydrogen peroxide. This suggests that the cholesterol ester conversion is the rate-limiting step. Interestingly, the hydrogen peroxide concentration reaches a maximum value toward the middle of the membrane and decreases on both ends, possibly due to higher mean cholesterol ester and cholesterol concentration in the middle region.

The current density plot shows that the current density reaches steady state at time 2T (17 seconds) correlating with the steady concentration of cholesterol ester and cholesterol concentration at this time point (**Figure 6**). Due to low cholesterol ester conversion, the cholesterol concentration at steady state is low, 0.0028 (0.5 μM). The current density plot also shows that the biosensor signal reaches a steady current at 2T (17 seconds) with a current magnitude of 0.0032 (0.64 μA/(cm2 )).

**Figure 3.** The concentration profile of cholesterol ester. The concentration profile shows that only 25% of the cholesterol ester has been converted to cholesterol at steady-state time and maintains at a constant concentration.

differential equations is casted into matrix form and solved using the MATLAB function PDEPE by assuming slab geometry. All units are non-dimensionalized. Note that *X* = 1 is the top of the membrane and *X* = 0 is the side of sensing electrode. For the simplified model, the

The solution of the system reaction-diffusion equation for cholesterol analysis is shown in **Figures 3**–**5**. The concentration profile shows that only 25% of the cholesterol ester has been converted to cholesterol at steady-state time and maintains at a constant concentration. On the other hand, the free cholesterol achieves 100% conversion to hydrogen peroxide. This suggests that the cholesterol ester conversion is the rate-limiting step. Interestingly, the hydrogen peroxide concentration reaches a maximum value toward the middle of the membrane and decreases on both ends, possibly due to higher mean cholesterol ester and cholesterol concen-

The current density plot shows that the current density reaches steady state at time 2T (17 seconds) correlating with the steady concentration of cholesterol ester and cholesterol concentration at this time point (**Figure 6**). Due to low cholesterol ester conversion, the cholesterol concentration at steady state is low, 0.0028 (0.5 μM). The current density plot also shows that the biosensor signal reaches a steady current at 2T (17 seconds) with a current magnitude of

**Figure 3.** The concentration profile of cholesterol ester. The concentration profile shows that only 25% of the cholesterol

ester has been converted to cholesterol at steady-state time and maintains at a constant concentration.

Thiele modulus for both CE and CO is 1.

128 Computer-aided Technologies - Applications in Engineering and Medicine

)).

tration in the middle region.

0.0032 (0.64 μA/(cm2

**Figure 4.** The concentration profile of free cholesterol. The free cholesterol achieves 100% conversion to hydrogen peroxide, suggesting that the cholesterol ester conversion is the rate-limiting step.

**Figure 5.** The concentration profile of the reaction product hydrogen peroxide. The hydrogen peroxide concentration reaches a maximum value toward the middle of the membrane and decreases on both ends, possibly due to higher mean cholesterol ester and cholesterol concentration in the middle region.

**Figure 6.** Plot of current density as a function of time derived from hydrogen peroxide reaction profile (**Figure 5**) at *X* = 0.

#### **3.4. Biosensor performance characterization**

#### *3.4.1. Substrate-current activity relationship*

The current-substrate calibration curve is constructed by simulating the steady-state current at 2T and varying the concentration of cholesterol ester from 1 to 10 S. The plot shows that the biosensor has high sensitivity and is linear at low substrate concentration (<1 S unit) and the current density continues to increase and gradually reaches a plateau at high substrate concentration (> 1 S unit). The maximum current density is 0.064 (12.8 μA/(cm2 )) (**Figure 7**).

**Figure 7.** Characterization of biosensor performance using substrate-current activity relationship. (A) The maximum current density is predicted at different substrate concentrations. Note that the substrate-current activity relationship obeys the Michaelis-Menten kinetic equation. (B) Inverse linear plot of the inverse current and substrate showed high linearity (*R*<sup>2</sup> = 0.987).

To quantitatively model how current density varies with substrate concentration, a mathematical model was constructed. Here, we observed that substrate-current activity relationship obeys the Michaelis-Menten kinetic equation. Specifically,

$$J = \frac{a \begin{bmatrix} S \end{bmatrix}}{b + \begin{bmatrix} S \end{bmatrix}} \tag{22}$$

*J* is the dimensionless current density, [*S*] is the cholesterol ester substrate concentration, and *a* and *b* are undetermined variables. At low substrate concentration ([*S*] ≪1), the equation can be approximated by a linear equation with the slope proportional to the ratio of *a* and *b*:

$$J = \frac{a}{b} \begin{bmatrix} S \end{bmatrix} \tag{23}$$

At high substrate concentration, the current density is constant, corresponding to maximum current density:

$$J = J\_{\text{max}} = a \tag{24}$$

The two variables can be solved by the inverse linear plot:

$$\frac{1}{J} = \frac{b + \begin{bmatrix} S \end{bmatrix}}{a \begin{bmatrix} S \end{bmatrix}} = \frac{b}{a} \begin{bmatrix} S \end{bmatrix} + \frac{1}{a} \tag{25}$$

By curve fitting, we determine *a* = 0.0066 and *b* = 0.965. Thus, we have recovered the substratecurrent activity relationship:

$$J = \frac{0.0066 \text{[S]}}{0.96 \text{S} + \text{[S]}} \tag{26}$$

#### *3.4.2. Sensitivity*

**Figure 6.** Plot of current density as a function of time derived from hydrogen peroxide reaction profile (**Figure 5**) at

The current-substrate calibration curve is constructed by simulating the steady-state current at 2T and varying the concentration of cholesterol ester from 1 to 10 S. The plot shows that the biosensor has high sensitivity and is linear at low substrate concentration (<1 S unit) and the current density continues to increase and gradually reaches a plateau at high substrate

**Figure 7.** Characterization of biosensor performance using substrate-current activity relationship. (A) The maximum current density is predicted at different substrate concentrations. Note that the substrate-current activity relationship obeys the Michaelis-Menten kinetic equation. (B) Inverse linear plot of the inverse current and substrate showed high

)) (**Figure 7**).

concentration (> 1 S unit). The maximum current density is 0.064 (12.8 μA/(cm2

*X* = 0.

linearity (*R*<sup>2</sup>

= 0.987).

**3.4. Biosensor performance characterization**

130 Computer-aided Technologies - Applications in Engineering and Medicine

*3.4.1. Substrate-current activity relationship*

The sensitivity of a biosensor signal can be defined as the unit increases in current density per unit increase in substrate concentration at the steady-state condition. High sensitivity is often required to differentiate between normal and pathological readings for a biosensor measurement. From the substrate-current activity relationship, the sensitivity can be expressed as the derivative of current density *J* with respect to [*S*]. To determine the optimal sensitivity, we calculate the second derivatives of the sensitivity expression at small *S* concentrations:

$$\frac{d^2J}{dS^2} = \frac{2a\left[S\right]}{\left(b + \left[S\right]\right)^3} - \frac{2a}{\left(b + \left[S\right]\right)^2} \equiv \frac{2a}{b + \left[S\right]^2} \left(\frac{\left[S\right]}{b} - 1\right) \tag{27}$$

Thus, we find that the substrate concentration [*S*] at *b* = 0.965 yields the greatest sensitivity based on the current density output.

#### *3.4.3. Measurement time*

The measurement time is another important element in biosensor design. For a point-of-care biosensor, it determines how fast the patient can retrieve the test results after the serum has been applied to the test strip. For an amperometric biosensor, the measurement time is determined by how fast the current achieves steady state. In our example, the measurement time is approximated at 2T or 17 seconds.

#### *3.4.4. Measurement range*

The measurement range is defined as the substrate concentration range where the substratecurrent activity relationship is linear. A broad measurement range is preferable to be able to differentiate normal from pathological condition as well as pathological severity. As the calibration curve often plateaus out at higher substrate concentration, it becomes harder to provide an accurate reading at this region. In our example, the linear region is at [*S*] < 2 (38 μM).

### **3.5. Biosensor optimization**

Since the behavior of the reaction-diffusion equation is controlled by two Thiele moduli ϕ1 and ϕ2 of CE and CO, it is possible to optimize biosensor response using these two parameters. From Eqs. (11) and (13), the conversion rate of cholesterol ester and the production rate of hydrogen peroxide are proportional to ϕ12 and ϕ22 , respectively. From the computational simulation, we observed that Thiele modulus 1 = 1 at *S* = 1 only achieved a 25% cholesterol ester conversion rate. By increasing the Thiele moduli 1–5, the system was able to achieve full cholesterol ester conversion at 1T (8.5 seconds). Similarly, by increasing the Thiele moduli 2– 5, the steady current attained a maximum value 0.05 (10 μA/(cm2 )). The results showed a 25 fold increase in sensitivity from the original value of 0.007 to 0.05, and the maximum current increased from 0.006 (1.2 μA/(cm2 )) to 0.55 (0.1 mA/(cm2 )) (**Table 1**).

**Table 1.** Biosensor optimization using the Thiele modulus.

#### **3.6. Biosensor model validation**

[ ]

132 Computer-aided Technologies - Applications in Engineering and Medicine

based on the current density output.

time is approximated at 2T or 17 seconds.

hydrogen peroxide are proportional to ϕ12

**Table 1.** Biosensor optimization using the Thiele modulus.

increased from 0.006 (1.2 μA/(cm2

5, the steady current attained a maximum value 0.05 (10 μA/(cm2

1 1 0.0068 1 1 5 5 0.053 1 1 8 16 0.05 20 1

*3.4.3. Measurement time*

*3.4.4. Measurement range*

**3.5. Biosensor optimization**

μM).

( [ ]) ( [ ]) [ ] [ ] <sup>2</sup> 2 3 22 <sup>2</sup> 2 2 <sup>1</sup> *d J a S a a <sup>S</sup> dS bS bS b S b*

=- @ - ç ÷

Thus, we find that the substrate concentration [*S*] at *b* = 0.965 yields the greatest sensitivity

The measurement time is another important element in biosensor design. For a point-of-care biosensor, it determines how fast the patient can retrieve the test results after the serum has been applied to the test strip. For an amperometric biosensor, the measurement time is determined by how fast the current achieves steady state. In our example, the measurement

The measurement range is defined as the substrate concentration range where the substratecurrent activity relationship is linear. A broad measurement range is preferable to be able to differentiate normal from pathological condition as well as pathological severity. As the calibration curve often plateaus out at higher substrate concentration, it becomes harder to provide an accurate reading at this region. In our example, the linear region is at [*S*] < 2 (38

Since the behavior of the reaction-diffusion equation is controlled by two Thiele moduli ϕ1 and ϕ2 of CE and CO, it is possible to optimize biosensor response using these two parameters. From Eqs. (11) and (13), the conversion rate of cholesterol ester and the production rate of

simulation, we observed that Thiele modulus 1 = 1 at *S* = 1 only achieved a 25% cholesterol ester conversion rate. By increasing the Thiele moduli 1–5, the system was able to achieve full cholesterol ester conversion at 1T (8.5 seconds). Similarly, by increasing the Thiele moduli 2–

fold increase in sensitivity from the original value of 0.007 to 0.05, and the maximum current

)) to 0.55 (0.1 mA/(cm2

**ϕ1 ϕ2 Sensitivity (***dJ***/***dS***) Measurement range (***S***) Measurement time (***T***)**

and ϕ22

, respectively. From the computational

)) (**Table 1**).

)). The results showed a 25-

æ ö

+ + <sup>+</sup> è ø (27)

The biosensor simulation result was compared with real-time experimental data [12]. The biosensor system consisted of a planar gold electrode modified with immobilized enzyme peroxidase, an insoluble mediator, layers of acetate cellulose, and the enzyme cholesterol oxidase in a buffer solution. Furthermore, all the cholesterol in the ester form was pretreated with cholesterol esterase to convert it to free cholesterol. The experiment showed that the linear range is up to 40 μM. We applied the optimization techniques to completely convert and linearize the cholesterol ester calibration up to 2S or 38 μM using a steady-state current. The optimized Thiele moduli were found to be Thiele modulus 1 = 4 and Thiele modulus 2 = 5, respectively. The membrane thickness was chosen at 0.1 um to match the maximum current in the experiment 70 nA and assumed that the reaction area diameter was 5mm. The simulation results showed excellent agreement with the experimental data. The calibration curve showed that the current response was linear up to 50 μM with a sensitivity of 0.5 nA/μM, which was close to the experiment values of the maximum range 40μM and sensitivity 0.6 nA/μM (**Figure 8**).

**Figure 8.** Biosensor model validation. The biosensor current response of cholesterol input predicted by the enzymatic diffusion equation is compared with that derived from the experiment. (A) Biosensor current response by mathematical simulation. Simulation parameters: l = 0.1, cholesterol esterase = 0.029 mg, cholesterol oxidase = 0.21mg, and the Biot number = ∞. (B) Biosensor current response by experimentation. Planar gold electrode modified with immobilized enzyme peroxidase, an insoluble mediator, layers of acetate cellulose, and the enzyme cholesterol oxidase in a buffer solution.

## **4. Conclusion**

We have constructed a mathematical model to simulate biosensor kinetics for cholesterol determination. By monitoring the steady current from the current density plot, the measurement time can be determined, and the current-substrate calibration can be constructed at the steady-state time point. The maximum range is the maximum concentration of cholesterol ester

in the linear region. The characteristic substrate-current activity relationship obeys Michaelis-Menten kinetics, and the equation can be determined by the sensitivity and maximum current density. Through mathematical analysis of system parameters, the optimal enzyme composition is shown to be the minimum amount of enzymes that can achieve total conversion of cholesterol ester as well as produce the highest possible steady-state current. Finally, the model was validated by comparing the results with real-time biosensor data. The free parameters such as membrane thickness or reaction area can be scaled to match specific design criteria or scale-up operation. Several potential improvements for the biosensor model include incorporating a serum diffusion layer by reducing the Biot number to produce a more realistic simulation, generalizing the multienzyme formulation and simulating the biosensor response under different amplification methods. Finally, we anticipated that the computer-aided model could be applied to simulate other multienzyme systems such as glucose or triglyceride determination as well as any other point-of-care diagnostic application.

## **Acknowledgements**

We thank Dr. Elliot Landaw at UCLA for the helpful discussion of the materials.

## **Author details**

Yu-Chen Lo1\*, Ren Gui2 , Hiroshi Honda2 and Jorge Z. Torres3

\*Address all correspondence to<EMAIL_ADDRESS>
1 Department of Bioengineering, Stanford University, Stanford, California, United States

2 Department of Bioengineering, Northwestern Polytechnic University, Fremont, California, United States

3 Department of Chemistry and Biochemistry, University of California, Los Angeles, California, United States

## **References**

[1] Kavey RE, Daniels SR, Lauer RM, Atkins DL, Hayman LL, Taubert K, et al. American Heart Association guidelines for primary prevention of atherosclerotic cardiovascular disease beginning in childhood. Circulation. 2003;107(11):1562-6. PubMed PMID: 12654618.

[2] Romero MR, Ahumada F, Garay F, Baruzzi AM. Amperometric biosensor for direct blood lactate detection. Analytical Chemistry. 2010;82(13):5568-72. doi: 10.1021/ ac1004426. PubMed PMID: 20518470.

in the linear region. The characteristic substrate-current activity relationship obeys Michaelis-Menten kinetics, and the equation can be determined by the sensitivity and maximum current density. Through mathematical analysis of system parameters, the optimal enzyme composition is shown to be the minimum amount of enzymes that can achieve total conversion of cholesterol ester as well as produce the highest possible steady-state current. Finally, the model was validated by comparing the results with real-time biosensor data. The free parameters such as membrane thickness or reaction area can be scaled to match specific design criteria or scale-up operation. Several potential improvements for the biosensor model include incorporating a serum diffusion layer by reducing the Biot number to produce a more realistic simulation, generalizing the multienzyme formulation and simulating the biosensor response under different amplification methods. Finally, we anticipated that the computer-aided model could be applied to simulate other multienzyme systems such as glucose or triglyceride

determination as well as any other point-of-care diagnostic application.

134 Computer-aided Technologies - Applications in Engineering and Medicine

, Hiroshi Honda2

\*Address all correspondence to<EMAIL_ADDRESS>
We thank Dr. Elliot Landaw at UCLA for the helpful discussion of the materials.

and Jorge Z. Torres3

1 Department of Bioengineering, Stanford University, Stanford, California, United States

2 Department of Bioengineering, Northwestern Polytechnic University, Fremont, California,

3 Department of Chemistry and Biochemistry, University of California, Los Angeles, Califor-

[1] Kavey RE, Daniels SR, Lauer RM, Atkins DL, Hayman LL, Taubert K, et al. American Heart Association guidelines for primary prevention of atherosclerotic cardiovascular disease beginning in childhood. Circulation. 2003;107(11):1562-6. PubMed PMID:

**Acknowledgements**

**Author details**

United States

nia, United States

**References**

12654618.

Yu-Chen Lo1\*, Ren Gui2

#### **Computer-Aided Diagnosis in Neuroimaging Computer-Aided Diagnosis in Neuroimaging**

Francisco J. Martínez-Murcia, Francisco J. Martínez-Murcia,

Juan Manuel Górriz and Javier Ramírez Juan Manuel Górriz and Javier Ramírez

Additional information is available at the end of the chapter Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/64980

#### **Abstract**

This chapter is intended to provide an overview to the most used methods for computeraided diagnosis in neuroimaging and its application to neurodegenerative diseases. The fundamental preprocessing steps, and how they are applied to different image modalities, will be thoroughly presented. We introduce a number of widely used neuroimaging analysis algorithms, together with a wide overview on the recent advances in brain imaging processing. Finally, we provide a general conclusion on the state of the art in brain imaging processing and possible future developments.

**Keywords:** neuroimaging, VBM, feature extraction, CT, MRI, PET, SPECT, machine learning, classification

## **1. Introduction**

Neuroimaging has meant a major breakthrough for the diagnosis and treatment of neurodegenerative diseases. Not so long ago, biomedical signal processing was limited to filtering, modelling or spectral analysis, prior to visual inspection. In the past decades, a number of powerful mathematical and statistical tools have been developed and evolved together with an increasing development and use of neuroimaging. Structural modalities such as computed tomography (CT) or the widely known magnetic resonance imaging (MRI), and later functional imaging techniques such as positron emission tomography (PET) and single photon emission computed tomography (SPECT) provide unprecedented insight in the internals of the brain, allowing the study of the structural and functional changes that can be linked to neurodegenerative diseases. This means a huge amount of data where automatic tools can help to identify patterns, reduce noise and enhance our knowledge of the brain functioning.

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Computer-aided diagnosis (CAD) systems in neuroimaging include a variety of methods that range from preprocessing of the images (just after acquisition) to advanced machine-learning algorithms to identify disease-related patterns. Algorithms used in the reconstruction of medical imaging, such as the tomographic reconstruction (TR) or the filtered back-projection (FBP) lay outside the scope of this chapter, focused on the application of CAD systems to neuroimaging.

This chapter starts with an exposition of the preprocessing methods used in different neuroimaging modalities, including registration, normalization and segmentation. We provide references on the algorithms behind well-known pieces of software such as statistical parametric mapping (SPM) [1], FreeSurfer [2] or the FMRIB Software Library (FSL) [3]. Later, the most used computer-aided diagnosis systems in psychiatry, psychology and neurology are described. These include SPM [1] and voxel-based morphometry (VBM) [4], voxels as features (VAFs) [5] and how the computation of regions of interest (ROIs) work in semiquantitative analysis. In the next section, new advances in neuroimaging analysis are presented, starting with the basis of machine learning and classification, including support vector machines (SVMs) [5–8], but also logistic regression [9, 10] or classifier ensembles [11, 12]. Given the characteristics of neuroimaging data, where we study large, possibly correlated, data, the extraction of higher-level features is essential. Therefore, in the last section, we provide an introduction to commonly used image decomposition algorithms such as principal component analysis (PCA) [8, 13–18] and independent component analysis (ICA) [19–22]. Finally, other recent feature extraction algorithms including spatial and statistical methods such as texture analysis [23–31], morphological tools [31–33] or artificial neural networks [34–40] are presented.

## **2. Preprocessing**

Preprocessing of the neurological images is a fundamental step in CAD systems as it ensures that all the images, either structural or functional, are comparable. We consider a preprocessing step, an algorithm that, applied after the acquisition and reconstruction of the images–usually a machine-dependent procedure–is intended to produce directly comparable images that represent a certain magnitude. The number and type of procedures to follow in preprocessing differs from one modality to another, although normalization and smoothing are used throughout all of them (see **Figure 1**).

## **2.1. Spatial normalization or registration**

The anatomy of every subject's brain is slightly different in shape and size. In order to compare images of different subjects, we need to eliminate these particularities and transform the images so that the subsequent group analysis or comparison can be performed. To do so, the individual images are mapped from their individual subject space (current anatomy) to a reference space, a common anatomical reference that allows the comparison. This procedure is known as *spatial normalization* or *registration*.

There are a number of algorithms used in image registration, but the procedure usually involves the computation of a series of parameters to map the source images to a template that works as a common anatomical reference (see Section 2.1.2 for an overview on registration algorithms). The most widely used template is the Montreal Neurological Institute (MNI) template.

**Figure 1.** Example of the pipeline followed in structural MRI preprocessing, comprising spatial normalization, smoothing and segmentation.

#### *2.1.1. The MNI space and template*

Computer-aided diagnosis (CAD) systems in neuroimaging include a variety of methods that range from preprocessing of the images (just after acquisition) to advanced machine-learning algorithms to identify disease-related patterns. Algorithms used in the reconstruction of medical imaging, such as the tomographic reconstruction (TR) or the filtered back-projection (FBP) lay outside the scope of this chapter, focused on the application of CAD systems to

138 Computer-aided Technologies - Applications in Engineering and Medicine

This chapter starts with an exposition of the preprocessing methods used in different neuroimaging modalities, including registration, normalization and segmentation. We provide references on the algorithms behind well-known pieces of software such as statistical parametric mapping (SPM) [1], FreeSurfer [2] or the FMRIB Software Library (FSL) [3]. Later, the most used computer-aided diagnosis systems in psychiatry, psychology and neurology are described. These include SPM [1] and voxel-based morphometry (VBM) [4], voxels as features (VAFs) [5] and how the computation of regions of interest (ROIs) work in semiquantitative analysis. In the next section, new advances in neuroimaging analysis are presented, starting with the basis of machine learning and classification, including support vector machines (SVMs) [5–8], but also logistic regression [9, 10] or classifier ensembles [11, 12]. Given the characteristics of neuroimaging data, where we study large, possibly correlated, data, the extraction of higher-level features is essential. Therefore, in the last section, we provide an introduction to commonly used image decomposition algorithms such as principal component analysis (PCA) [8, 13–18] and independent component analysis (ICA) [19–22]. Finally, other recent feature extraction algorithms including spatial and statistical methods such as texture analysis [23–31], morphological tools [31–33] or artificial neural networks [34–40] are

Preprocessing of the neurological images is a fundamental step in CAD systems as it ensures that all the images, either structural or functional, are comparable. We consider a preprocessing step, an algorithm that, applied after the acquisition and reconstruction of the images–usually a machine-dependent procedure–is intended to produce directly comparable images that represent a certain magnitude. The number and type of procedures to follow in preprocessing differs from one modality to another, although normalization and smoothing are used

The anatomy of every subject's brain is slightly different in shape and size. In order to compare images of different subjects, we need to eliminate these particularities and transform the images so that the subsequent group analysis or comparison can be performed. To do so, the individual images are mapped from their individual subject space (current anatomy) to a reference space, a common anatomical reference that allows the comparison. This procedure

neuroimaging.

presented.

**2. Preprocessing**

throughout all of them (see **Figure 1**).

**2.1. Spatial normalization or registration**

is known as *spatial normalization* or *registration*.

The MNI space is the most widely used space for brain registration and was recently adopted by the International Consortium for Brain Mapping (ICBM) as its standard template. It defines a standard three-dimensional (3D) coordinate system (also known as 'atlas'), which is used to map the location of brain structures independently of the size and shape of each subject's brain.

The MNI space was intended to replace the Talairach space, a system based on a dissected and photographed brain for the Talairach and Tournoux atlas. In contrast to this, the MNI created a new template that was approximately matched to the Talairach brain but using a set of normal MRI scans. The current standard MNI template is the ICBM152 [41], which is the average of 152 normal MRI scans that have been matched to an older version of the MNI template using a nine-parameter affine transform.

## *2.1.2. Registration algorithms*

Algorithms used in registration can be categorized in *linear* transformations (being the *affine* transform the most complex) and *non-rigid* or elastic transformations. Affine transformations are applied as a matrix multiplication and include terms for translation, scale, rotation, shear, squeeze and so on. These lineal transformations are applied globally to the image and therefore do not account for local geometric differences. Most neuroimaging software includes some kind of affine registration, including FreeSurfer [2], FSL (via package FLIRT) [3] or Elastix.

The estimation of the parameters is performed via the optimization of a given cost function, the minimum-squared difference between the source image and the template being the most basic. Modern software include more refined functions, for example, Tukey's biweight function (in mri\_robust\_template of FreeSurfer) [11], or the mutual information (in FLIRT) [42], that operate under a high-complexity schema involving local and global multiresolution optimization. When working with images of the same modality, the preferred cost function is the minimum-squared difference between the source image and the template, whereas in the case of multimodal registration, the maximization of the mutual information is preferred.

Non-rigid transformations can apply local transformation to align the target image to the template. Many of these non-rigid transformations are applied as a local fine-tuning after a previous affine transformation, although some of them use higher-complexity models that do not need this previous step. Some non-rigid transformations include radial basis functions (RBFs) (thin plate or surface splines, multiquadrics and compactly supported transformations), physical continuum models and large deformation models (diffeomorphisms). Of these, the most popular are diffeomorphic transformations, which feature the estimation and application of a warp field, and they are used in SPM (default) [1], FreeSurfer [2], FSL (FNIRT) [3], Elastix or ANTs.

## *2.1.3. Co-registration*

Sometimes, we use different modalities, usually a functional and a structural image, for the same subject. It would be therefore very useful that these two (or even more) different images spatially match each other, so that any processing can be fit to any of them and applied to all. Since functional images have low resolution, the procedure for performing spatial normalization frequently involves co-registration to the structural image, which has a higher level of detail. Then the spatial normalization (or registration) parameters are estimated on the structural image and applied to all the co-registered maps. In the context of low-resolution functional imaging, affine co-registration to the structural image is preferred.

## **2.2. Smoothing**

Despite the spatial normalization applied to the images, differences between subjects are still a problem that can reduce the signal-to-noise ratio (SNR) of the neuroimaging data. This problem increases as the number of subjects increases. To increase the SNR, it is recommended to filter out the highest frequencies, that is, applying a *smoothing*. Smoothing removes the smallest scale changes between voxels, making the detection of large-scale changes easier.

On the other hand, smoothing the images lowers its resolution. Smoothing is usually applied using a 3D Gaussian kernel to the image, determined by its full-width at half maximum (FWHM) parameter, which is the diameter of the smoothing kernel at half of its height. The choice of the size of the kernel is therefore of fundamental importance, and it depends on the signal to be detected. A small kernel will make our further processing confound noise with activation signal, but a larger kernel can parse out significant signal. As a general rule, the size of the kernel should be smaller than the activation to be detected.

## **2.3. Functional MRI-specific steps**

The estimation of the parameters is performed via the optimization of a given cost function, the minimum-squared difference between the source image and the template being the most basic. Modern software include more refined functions, for example, Tukey's biweight function (in mri\_robust\_template of FreeSurfer) [11], or the mutual information (in FLIRT) [42], that operate under a high-complexity schema involving local and global multiresolution optimization. When working with images of the same modality, the preferred cost function is the minimum-squared difference between the source image and the template, whereas in the case

140 Computer-aided Technologies - Applications in Engineering and Medicine

of multimodal registration, the maximization of the mutual information is preferred.

or ANTs.

*2.1.3. Co-registration*

**2.2. Smoothing**

Non-rigid transformations can apply local transformation to align the target image to the template. Many of these non-rigid transformations are applied as a local fine-tuning after a previous affine transformation, although some of them use higher-complexity models that do not need this previous step. Some non-rigid transformations include radial basis functions (RBFs) (thin plate or surface splines, multiquadrics and compactly supported transformations), physical continuum models and large deformation models (diffeomorphisms). Of these, the most popular are diffeomorphic transformations, which feature the estimation and application of a warp field, and they are used in SPM (default) [1], FreeSurfer [2], FSL (FNIRT) [3], Elastix

Sometimes, we use different modalities, usually a functional and a structural image, for the same subject. It would be therefore very useful that these two (or even more) different images spatially match each other, so that any processing can be fit to any of them and applied to all. Since functional images have low resolution, the procedure for performing spatial normalization frequently involves co-registration to the structural image, which has a higher level of detail. Then the spatial normalization (or registration) parameters are estimated on the structural image and applied to all the co-registered maps. In the context of low-resolution

Despite the spatial normalization applied to the images, differences between subjects are still a problem that can reduce the signal-to-noise ratio (SNR) of the neuroimaging data. This problem increases as the number of subjects increases. To increase the SNR, it is recommended to filter out the highest frequencies, that is, applying a *smoothing*. Smoothing removes the smallest scale changes between voxels, making the detection of large-scale changes easier.

On the other hand, smoothing the images lowers its resolution. Smoothing is usually applied using a 3D Gaussian kernel to the image, determined by its full-width at half maximum (FWHM) parameter, which is the diameter of the smoothing kernel at half of its height. The choice of the size of the kernel is therefore of fundamental importance, and it depends on the signal to be detected. A small kernel will make our further processing confound noise with activation signal, but a larger kernel can parse out significant signal. As a general rule, the size

functional imaging, affine co-registration to the structural image is preferred.

of the kernel should be smaller than the activation to be detected.

The acquisition of functional MRI (fMRI) involves a more complex preprocessing of the images, given its dynamic features. fMRI studies acquire a long sequence of low-resolution MRI images that contain a magnitude known as blood-oxygen-level-dependent (BOLD) contrast. This sequence of three-dimensional volumes is combined into a four-dimensional (4D) volume that conceptually works similar to a video. In this context, the outcome can be much more affected by the subject motion. Therefore, procedures such as slice-timing correction and motion correction are mandatory in fMRI.

## *2.3.1. Slice-timing correction*

In fMRI, the scanner acquires each slice within a single brain volume at different times. Different methods are used to acquire the slices: descending order (top to down), ascending order (bottom to up) and interleaved (the slices are acquired in a certain sequence). The time interval between one slice and another is usually small, but after acquiring a whole brain volume, there might be a difference of several seconds between the first and the last acquisition.

To compensate for the time differences between slice acquisition within a single volume, a *slice-timing correction* is performed by temporally interpolating the slices so that the resulting volume is approximately equivalent to acquiring the whole brain at the same time point. These time differences are especially important in event-driven experimental design, where timing is very relevant. Linear, quadratic or spline functions are used to interpolate the slices.

## *2.3.2. Motion correction*

*Motion correction*, also known as *realignment*, accounts for the undesirable effect of head movement during the acquisition of the data. It is almost impossible for a subject to lie perfectly still during the acquisition of an fMRI sequence, and even the small movements can lead to variation in the BOLD values that, if uncorrected, can lead to activation of false positives.

The correction of this problem usually uses a rigid-body transformation similar to those used in registration. In this case, a model characterized by six parameters that account for translation and rotation is frequently used. The parameter estimation is performed by minimizing a cost function, such as correlation or mutual information, between the volumes and a reference time volume, usually the mean image of all time points.

Sometimes, the movement of the head is so fast that motion correction cannot correct its effects. In that case, the most used approach is to eliminate the images acquired during that fast movement, using an artefact detection algorithm that identifies large variations between images at adjacent time points.

#### **2.4. Intensity normalization**

Most functional neuroimaging modalities, in contrast to unitless structural MRI images, are the representation of the distribution of a certain contrast over the brain. There exist a larger number of sources of variability that can affect the final values: contrast uptake, radiotracer decay time, metabolism, and so on. In order to establish comparisons between subjects, an *intensity normalization* procedure is required.

Intensity normalization methods are to be linear in nature, since it is essential to maintain the intensity ratio between brain regions, acting on the whole brain. In its simplest form, it consists of a division by a constant. This parameter is often estimated [6, 7] as the average value of the 95th bin of the histogram of the image, that is, the average of the 5% higher-intensity values, in what is known as the *normalization to the maximum*. Another approach, called *integral normalization*, estimates this parameter as the sum of all values of the image. A more complex approach requires *a priori* knowledge of the distribution of intensity in a normal subject. This is designed so that the whole image is divided by the binding potential (BP) [43], a specific-tonon-specific ratio between the intensities in areas where the tracer should be concentred and the non-specific areas.

Then, we have general linear transformations, defined as *Y* = *aX* + *b*. These procedures use estimates of the probability density function (PDF) of the source images and then estimate the parameters *a* and *b*, which transform their original PDF to an expected range. Methods to estimate the PDF parameters range from the simplest, non-parametric histogram [44] to more advanced estimates such as an analysis of covariance (ANCOVA) approach used in SPM [44, 45], or parametric estimates involving the Gaussian or the more general alpha-stable distribution [46], which has been recently tested with great success.

Structural modalities also suffer from some sources of intensity variability, for example, magnetic field inhomogeneity, noise, evolution of the scanners, and so on. Field inhomogeneity causes distortions in both geometry and intensity of the MR images [47], usually addressed via increasing the strength of the gradient magnetic field or preprocessing. Intensity variability is especially noticeable in multicentre-imaging studies, where images should share certain characteristics. To improve the homogeneity of a set of structural images acquired at different locations, the use of quantitative MRI images has been recently proposed [48]. In contrast to typical unitless *T*1-weighted images, quantitative imaging can provide neuroimaging biomarkers for myelination, water and iron levels that are absolute measures comparable across imaging sites and time points.

#### **2.5. Segmentation**

Segmentation, mostly of structural MRI images, involves a series of algorithms aimed at constructing maps of the distribution of different tissues. The general approach is to separate the image in three different maps containing grey matter (GM), white matter (WM) and cerebrospinal fluid (CSF), although some software can also output data for bone, soft tissue or very detailed functional regions and subregions [49–51]. The procedure is applied after all aforementioned steps, including field inhomogeneity correction, which is essential for a correct segmentation. Here, we provide insight on the most used algorithms for segmentation

#### *2.5.1. Statistical parametric mapping*

decay time, metabolism, and so on. In order to establish comparisons between subjects, an

Intensity normalization methods are to be linear in nature, since it is essential to maintain the intensity ratio between brain regions, acting on the whole brain. In its simplest form, it consists of a division by a constant. This parameter is often estimated [6, 7] as the average value of the 95th bin of the histogram of the image, that is, the average of the 5% higher-intensity values, in what is known as the *normalization to the maximum*. Another approach, called *integral normalization*, estimates this parameter as the sum of all values of the image. A more complex approach requires *a priori* knowledge of the distribution of intensity in a normal subject. This is designed so that the whole image is divided by the binding potential (BP) [43], a specific-tonon-specific ratio between the intensities in areas where the tracer should be concentred and

Then, we have general linear transformations, defined as *Y* = *aX* + *b*. These procedures use estimates of the probability density function (PDF) of the source images and then estimate the parameters *a* and *b*, which transform their original PDF to an expected range. Methods to estimate the PDF parameters range from the simplest, non-parametric histogram [44] to more advanced estimates such as an analysis of covariance (ANCOVA) approach used in SPM [44, 45], or parametric estimates involving the Gaussian or the more general alpha-stable distri-