Guanylate-binding proteins consist of an N-terminal α, β globular large GTPase domain and a α-helical finger-like C-terminal regulatory domain. The domains are connected by a short intermediate region consisting of one α-helix and a short two-stranded β-sheet (12, 13). A GTPase-domain binds guanine nucleotides with low affinities. This induces nucleotide dependent GBP multimerization and hydrolysis of GTP via GDP to GMP [reviewed in Ref. (3)]. Human GBP1, GBP2, and GBP5 and murine Gbp2b/Gbp1, Gbp2, and Gbp5 have at the C-terminus a CaaX sequence (C—cysteine, aa two amino acids, X—terminal amino acid), which directs isoprenylation—the addition of lipid moiety to the protein, which targets proteins to intracellular membranes and facilitates protein-protein interaction (4). Recruitment of proteins to parasitophorous vacuoles harboring pathogens can lead to restriction of pathogen proliferation (14).

GBPs are involved in regulation of inflammasomes—a highmolecular-weight complexes present in the cytosol of stimulated immune cells that mediate the activation of inflammatory caspases resulting in pathogen clearance and/or death of infected cell [reviewed in Ref. (1, 3, 15)]. Gbps can also attack parasites directly via supramolecular complexes (16) and interfere with virus replication (17) or virion assembly (18). Type of effective defense depends on pathogen involved.

A wide range of studies revealed an important role of GBPs in response to different infections including viral (17–20), bacterial (21–24), and protozoan pathogens (14, 16, 25), both vacuolar (14, 16, 21, 24, 25) and cytosolic (17–20).

For example, in human GBP1 influences resistance to vesicular stomatitis virus (19), encephalomyelocarditis virus (19), influenza A viruses (17), and *Chlamydia trachomatis* (22), GBP3 reduces virus titers of influenza A viruses (17) and GBP5 prevents processing and incorporation of the viral glycoprotein Env of HIV-1 (18).

Murine Gbp2b/Gbp1 plays role in defense against *Listeria monocytogenes* and *Mycobacterium bovis* BCG (23), Gbp2 inhibits replication of vesicular stomatitis virus and encephalomyelocarditis virus (20), *Toxoplasma gondii* (14), and *Salmonella typhimurium* (24), and Gbp5 protects against *S. typhimurium* (21) and *M. bovis* BCG (23). Moreover, several Gbps can cooperate for more effective defense. Gene specific-silencing using siRNA established that murine Gbp2b/Gbp1, Gbp5, Gbp7, and Gbp6/10 protect against *M. bovis* BCG and *L. monocytogenes.* A combination of siRNAs exacerbated the loss of function, which indicated that protective Gbps functioned cooperatively (23). Similarly, mutual molecular interactions of murine Gbp2b/Gbp1, Gbp2, Gbp3, Gbp5, and Gbp6 protected against *T. gondii* (16).

*Leishmania* is an obligatory intracellular mammalian pathogen that enters skin by the bite of female phlebotomine sand flies and infects so-called professional phagocytes (neutrophils, monocytes, and macrophages), as well as dendritic cells and fibroblasts. The major host cell is the macrophage where parasites reside inside parasitophorous vacuole, multiply, eventually rupturing the cell and spread to uninfected cells. Infected cells can spread to lymph nodes, spleen, liver, bone marrow, and sometimes lungs [reviewed in Ref. (26)]. The infection can remain asymptomatic or result in one of three main clinical syndromes: the cutaneous form of the disease in dermis, which can be localized or diffuse; mucocutaneous leishmaniasis in the mucosa and the visceral leishmaniasis characterized by splenomegaly and hepatomegaly that results from the metastatic spread of infection to the spleen and liver (27, 28). Manifestations of the disease depend on the infecting species, environmental and social factors, and the genotype of the mammalian host [reviewed in Ref. (26)].

There is very little known about a possible role of GBPs in *Leishmania* infection. Analysis of global gene expression of bone marrow derived macrophages from BALB/c mouse demonstrated upregulation of expression of *Gbp2b/Gbp1*, *Gbp2*, *Gbp3*, *Gbp6*, and *Gbp7* after 24 hours of infection with *Leishmania major* promastigotes (29). Dendritic cells generated from blood of healthy human donors exhibited increased expression of *GBP1* and *GBP2* 16 hours after infection by *L. major* promastigotes, whereas dendritic cells infected by *Leishmania donovani* had increased expression of *GBP1* (30). Comparison of global gene expression in skin lesions of *Leishmania braziliensis*-infected patients with skin of normal skin biopsies revealed upregulation of *GBP5* mRNA (31).

For our analysis, we selected two murine *Gbps* with the C-terminal CaaX sequence enabling targeting proteins to parasitophorous membranes (4). We studied expression of *Gbp2b/ Gbp1* and *Gbp5 in vivo* before and 8 weeks after *L. major* infection in 10 mouse strains from two genetically distant but internally related groups: CcS/Dem (BALB/c, STS, CcS-5, CcS-16, CcS-20) and OcB/Dem (O20, C57BL/10 (B10), C57BL/10-*H2pz*

**Abbreviations:** Gbp, guanylate-binding protein (murine); GBP, guanylate-binding protein (human); *Gbp2b*/*Gbp1*, murine gene coding this guanylate-binding protein was originally named *Gbp1* and later renamed *Gbp2b*; RCS, recombinant congenic strains; CcS—BALB/c-c-STS, series of recombinant congenic containing random 12.5% of genome of the donor strain STS/A (STS) on 87.5% genome of the background strain BALB/cHeA (BALB/c); OcB—O20-c-C57BL/10-*H*-*2pz* (B10.O20/ Dem), series of recombinant congenic containing random 12.5% (or 6.25% or less) of genome of the donor strain B10.O20/Dem (B10.O20) on the background strain O20/A (O20).

(B10.O20), OcB-9, OcB-43). Each CcS/Dem strain contains a different, random set of approximately 12.5% genes of the donor strain STS and approximately 87.5% genes of the background strain BALB/c (32, 33). OcB/Dem strains were derived from strains B10.O20 and O20. Strains OcB-43 and OcB-9 contain different 4 and 12.5% of B10 genome on O20 background, respectively; strain B10.O20 contains 4% of O20 genome on B10 background (32, 33). The limited and defined genetic differences between strains in each group (33) make it possible to identify the differences in *Gbp* expression that are controlled by genes outside the *Gbp* coding gene-complex on chromosome 3. Incidence and size of skin lesions indicate that BALB/c and CcS-16 are highly susceptible and B10.O20 is susceptible to *L. major*; whereas CcS-20 is intermediate and STS, CcS-5, O20, B10, OcB-9, and OcB-43 are resistant to this parasite (34) (this study).

We found that the levels of *Gbp2b/Gbp1* and *Gbp5* mRNAs are influenced by *L. major* infection and by genome of the host. The infection caused a large increase of *Gbp2b/Gbp1* and *Gbp5* expression, but *Gbp*s levels in both uninfected and infected mice differ among mouse strains indicating influence of genetic factors. These genetic influences are different in uninfected and infected mice and in some strains there is a clear evidence for a regulation by genes other than the *Gbp* genes (*trans*-regulation). We also show that Gbp2b/Gbp1 protein and *L. major* parasites co-localize in resistant strains STS, CcS-5, and O20 and in the intermediate strain CcS-20 but not in the highly susceptible strain BALB/c.

#### MATERIALS AND METHODS

#### Mice

#### mRNA Expression Experiments

A total of 275 (152 infected and 123 uninfected) female mice of strains BALB/c (33 infected and 22 uninfected), STS (20 infected and 13 uninfected), CcS-5 (11 infected and 10 uninfected), CcS-16 (10 infected and 11 uninfected), CcS-20 (12 infected and 12 uninfected), O20/A (O20) (12 infected and 12 uninfected), C57BL/10Sn (B10) (17 infected and 10 uninfected), B10.O20/ R164/Dem (B10.O20) (17 infected and 12 uninfected), OcB-9 (7 infected and 7 uninfected), and OcB-43 (13 infected and 14 uninfected) were tested in 15 independent experiments. The age of mice was 8–23 weeks (mean = 11.9 weeks, median = 11 weeks) at the time of infection (start of experiment in control mice). A total of 81 infected mice of strains BALB/c (*n*= 9), STS (*n*= 10), CcS-5 (*n*= 11), O20 (*n*= 12), B10 (*n*= 16), B10.O20 (*n*= 16), and OcB-9 (*n* = 7) from these experiments were used for estimation of parasite load in skin and/or spleen. 40 infected female mice of the strains BALB/c (*n* = 5), STS (*n* = 10), CcS-5 (*n* = 4), CcS-16 (*n* = 12), and CcS-20 (*n* = 9) from additional four experiments were also used for the estimation of parasite load in skin and/or spleen.

#### Immunohistochemistry Experiments

97 (48 infected and 49 uninfected) female mice of strain BALB/c (9 infected and 9 uninfected), STS (9 infected and 9 uninfected), CcS-5 (8 infected and 8 uninfected), CcS-20 (11 infected and 11 uninfected), and O20 (11 infected and 12 uninfected) were tested in two independent experiments. The age of mice was 8–18 weeks (mean 13 weeks, median 14 weeks) at the time of infection.

#### Ethics Statement

All experimental protocols utilized in this study comply with the Czech Government Requirements under the Policy of Animal Protection Law (No. 246/1992) and with the regulations of the Ministry of Agriculture of the Czech Republic (No. 207/2004), which are in agreement with all relevant European Union guidelines for work with animals and were approved by the Institutional Animal Care Committee of the Institute of Molecular Genetics AS CR and by Departmental Expert Committee for the Approval of Projects of Experiments on Animals of the Academy of Sciences of the Czech Republic (permissions Nr. 190/2010; 232/2012).

#### Parasite

*Leishmania major* LV 561 (MHOM/IL/67/LRCL 137 JERICHO II) was maintained in rump lesions of BALB/c females. Amastigotes were transformed to promastigotes using SNB-9 (35). 107 promastigotes from the passage two cultivated for 6 days were inoculated in 50 µl sterile saline s.c. into mouse rump (36). Control uninfected mice were injected by 50 µl sterile saline.

#### Disease Phenotype

The size of the skin lesions was measured every week using the Profi LCD Electronic Digital Caliper Messschieber Schieblehre Messer (Shenzhen Xtension Technology Co., Ltd. Guangdong, China), which has accuracy 0.02 mm. The mice were killed 8 weeks after inoculation. Skin, spleen, liver, and inguinal lymph nodes were collected for later analysis. Preparation of skin samples: approximately 3 mm of border skin surrounding lesion was taken. Hair was removed with scissors. A half of each skin sample was snap-frozen in liquid nitrogen for further RNA and DNA isolations. Another half was fixed in 4% formaldehyde for further paraffin embedding and immunohistochemical analysis. Samples from uninfected mice were obtained from the same rump area and used as a negative control.

## Quantification of Parasite Load by PCR-ELISA

Parasite load was measured in DNA from frozen skin and spleen samples using PCR-ELISA according to the previously published protocol (37). Briefly, total DNA was isolated using a TRI reagent (Molecular Research Center, Cincinnati, OH, USA) standard procedure (https://www.mrcgene.com/wp-content/ uploads/2014/06/TRI-LSMarch2017.pdf) or a modified proteinase K procedure (37). For PCR, two primers (digoxigenin-labeled F 5′-ATT TTA CAC CAA CCC CCA GTT-3′ and biotin-labeled R 5′-GTG GGG GAG GGG CGT TCT-3′; VBC Genomics Biosciences Research, Austria) were used for amplification of the 120-bp conservative region of the kinetoplast minicircle of *Leishmania* parasite, and 50 ng of extracted DNA was used per each PCR reaction. For a positive control, 20 ng of *L. major* DNA per reaction was amplified as a highest concentration of standard. A 33-cycle (expression experiments) or 26-cycle (immunostaining experiments) PCR reaction was used for quantification of parasites. Under these conditions, the amount of PCR product is linearly proportional to number of parasites (37). PCR product was measured by the modified ELISA (Pharmingen, San Diego, CA, USA). Concentration of *Leishmania* DNA was determined using the ELISA Reader Tecan and the curve fitter program KIM-E (Schoeller Pharma, Prague, Czech Republic) with least squares-based linear regression analysis.

#### RNA Isolation

Mouse spleens, skins, liver, and inguinal lymph nodes were snapped frozen by liquid nitrogen immediately after dissection and stored at −80°C until total RNA extraction. At the time of RNA isolation tissue were homogenized in TRI Reagent (Sigma-Aldrich, Inc., St. Louis, MO, USA) using Polytron PT 2100 homogenizer (Kinematica Ag, Luzern, Switzerland) and immediately followed by total RNA isolation according to the manufacturer's protocol. RNA concentration was measured with Nanodrop (NanoDrop Technologies, LLC, Wilmington, DL), and quality of RNA was estimated also using Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA). The isolated RNA was stored at −80°C.

#### Real-time PCR

One microgram of total RNA was diluted in 8 µl of sterile RNaseand DNase-free water, was treated with 1 µl DNase I (1 U/μl) and 1 µl DNase I reaction buffer (Promega Corporation, Madison, WI, USA), and used for subsequent reverse transcription. Singlestrand cDNA was prepared from total RNA using Promega firststrand synthesis method. DNase I-treated RNA was incubated for 10 minutes at 65°C, then cooled quickly on ice for 5 minutes, and then treated with 1 µl DNase I stop solution (Promega Corporation, Madison, WI, USA). For the next step, a mixture containing 1 µl of random hexamers primers (100 ng/1 µl) (Invitrogen, Carlsbad, CA, USA), 5 µl (50 ng/µl) of dNTP mix (Invitrogen, Carlsbad, CA, USA), 5 µl of the reaction buffer (Promega Corporation, Madison, WI, USA), 2.5 µl of RNase/ DNase-free water (Invitrogen, Carlsbad, CA, USA), and 0.5 µl of M-MLV Reverse Transcriptase RNAase H Minus Point Mutant (100 U/1 μl) (Promega Corporation, Madison, WI, USA) was added and followed by 60 minutes at 37°C. Single-strand cDNA was kept at −80°C until RT-PCR analysis. Real-time PCR was performed using a BioRad iQ iCycler Detection System (Bio-Rad Laboratories, Inc., Hercules, CA, USA). Primers were designed using Roche Universal ProbeLibrary, ProbeFinder version 2.45 for mouse (*Gbp2b*/*Gbp1*-F AAACCAGGAGGCTACTACCTTTTT, *Gbp2b/Gbp1*-R GTATTTTCTCAGCATCACTTCAGC; *Gbp5*-F TTCACCCAATCTAAGACCAAGAC, *Gbp5*-R AGCACCAG GCTTTCTAGACG; *Gapdh*-F AGAACATCATCCCTGCAT CC, *Gapdh*-R ACATTGGGGGTAGGAACAC). Reaction was performed in total volume of 25 µl, including 12.5 µl of 2× SYBR Green Supermix (Bio-Rad Laboratories, Inc., Hercules, CA, USA), 1 µl of each primer of *Gbp2b*/*Gbp1* and *Gbp5* genes (final concentration 6.6 µM), 7.5 µl of water (Invitrogen, Carlsbad, CA, USA), and 3 µl of the cDNA template. Six different samples from each experimental group were used, and all samples were tested in triplets. The average Ct values (cycle threshold) were used for quantification, and the relative amount of each mRNA was normalized to the housekeeping gene, *Gapdh* mRNA. Using the delta Ct value, relative expression was calculated [ratio (reference/ target) = 2 Ct (reference) − Ct (target)] × 10,000. All experiments included negative controls containing water instead of cDNA.

# Genotyping of *Gbp* Cluster in OcB Series

DNA was isolated from tails using a standard proteinase procedure. Strains O20, B10, B10.O20, OcB-9, and OcB-43 were genotyped using microsatellite markers D3Mit160 (size of B10 allele 137 bp, size of O20 allele 127 bp) and D3Mit17 (B10 allele 200 bp, O20 allele 174 bp) (Generi Biotech, Hradec Králové, Czech Republic): The DNA genotyping by PCR was performed as described elsewhere (38).

#### Immunostaining

After deparaffinization and rehydration, the 3 µm thick slices of skin tissue were 15 minutes heat induced in Tris-EDTA buffer (10 mM Tris, 1 mM EDTA, pH 8.5) for antigen retrieval. For fluorescent labeling of *Leishmania* parasite was used anti-*Leishmania* lipophosphoglycan mouse monoclonal antibody (cat. no. CLP003A; Cedarlane, Hornby, Canada) and TRITC-labeled IgM (115-025-020; Jackson ImmunoResearch, West Grove, PA) all diluted 1:500. Gbp2b/Gbp1 protein was stained with rabbit anti-Gbp1 Polyclonal antibody (PA5-23509; Thermo Fisher Scientific, Rockford, IL, USA) diluted 1:100 and anti-rabbit-AlexaFluor-647 (cat. no. 711-605-152; Jackson ImmunoResearch, West Grove, PA) diluted 1:500. Nuclei of the cells were stained with bisBenzimide H33258 (Sigma-Aldrich, St. Louis, MO, USA) 10 mg per 1 ml diluted 1:1,000. Images were captured with microscope Leica DM6000 objective HCX PL Apo 40×/0.75 PH2 and color camera Leica DFC490. Evaluation of images was done with Fiji ImageJ 1.51n. 10 fields (320.66 × 239.57 µM) from each mouse were analyzed.

#### Statistical Analysis

The differences among strains within each of the two groups in expression of *Gbp2b*/*Gbp1* and *Gbp5* and the differences between uninfected and infected mice were evaluated by Mann–Whitney test using the program Statistica for Windows 12.0 (StatSoft, Inc., Tulsa, OK, USA). The results were corrected for multiple testing by Bonferroni correction. The correction factor was 10× both for intragroup differences and differences between infected and uninfected mice of the same strain.

#### RESULTS

## Mouse Strains Differ in Expression of Both *Gbp2b*/*Gbp1* and *Gbp5* in Uninfected Mice

We observed strong genetic influence on mRNA levels of tested *Gbp*s. We have examined expression of *Gbp2b*/*Gbp1* (**Figure 1**) and *Gbp5* (**Figure 2**) in skin, inguinal lymph nodes, spleen, and liver of uninfected mice belonging to two genetically different series of strains—CcS/Dem (BALB/c, STS, CcS-5, CcS-16, CcS-20) and OcB/Dem (O20, B10, B10.O20, OcB-9, OcB-43). We have compared expression in parental strains BALB/c and STS with the

FIGURE 1 | Differences in expression of *Gbp2b*/*Gbp1* in organs of uninfected mice. Expression of *Gbp2b*/*Gbp1* in skin (A), lymph nodes (B), spleen (C), and liver (D) of uninfected female mice of strains BALB/c (*n* = 7 skin, 9 lymph nodes, 11 spleen, 9 liver), STS (6 skin, 9 lymph nodes, 8 spleen, 6 liver), CcS-5 (6 skin, 6 lymph nodes, 6 spleen, 6 liver), CcS-16 (6 skin, 6 lymph nodes, 6 spleen, 6 liver), CcS-20 (6 skin, 6 lymph nodes, 7 spleen, 6 liver), O20 (7 skin, 6 lymph nodes, 9 spleen, 6 liver), B10 (6 skin, 7 lymph nodes, 6 spleen, 6 liver), B10.O20 (5 skin, 7 lymph nodes, 6 spleen, 8 liver), OcB-9 (6 skin, 6 lymph nodes, 6 spleen, 5 liver), and OcB-43 (6 skin, 6 lymph nodes, 6 spleen, 7 liver) were compared. The data show the means ± SD. Only the differences between parental strains BALB/c and STS and strains of CcS/Dem series and parental strains O20 and B10 and strains of OcB/Dem series are shown. Nominal *P* values are shown.

strains of CcS/Dem series, and expression in parental strains O20 and B10 with the strains of OcB/Dem series (**Figures 1** and **2**), as well as expression of the strains within CcS/Dem and OcB/Dem series in skin (**Tables 1A,C**), lymph nodes (**Tables 2A,C**), spleen (**Tables 3A,C**), and liver (**Tables 4A,C**).

Expression of *Gbp2b*/*Gbp1* in background strain BALB/c and donor strain STS in skin (**Figure 1A**; **Table 1A**) does not differ, whereas strains CcS-5, CcS-16, and CcS-20, each carrying a different set of 12.5% genes of STS on BALB/c background, exhibit higher expression than either parent (**Figure 1A**; **Table 1A**). Expression of *Gbp2b*/*Gbp1* in parental strains of OcB/Dem series O20 and B10 in skin differed (*P* = 0.0047); strains B10.O20, OcB-9, and OcB-43 exceeded in *Gbp2b*/*Gbp1*

expression of the parental strain B10 but not O20 (**Figure 1A**; **Table 1A**).

We have observed differences in the expression of *Gbp2b*/*Gbp1* among mouse strains also in other tested organs (**Figures 1B–D**; **Tables 2A**, **3A**, and **4A**). Strains OcB-43 and OcB-9 differed in the expression of *Gbp2b*/*Gbp1* in lymph nodes (**Table 2A**), CcS-16 exhibited lower expression than STS in spleen (**Figure 2C**; **Table 3A**), B10 and B10.O20 exhibited lowest expression in spleen and differed from strains O20, OcB-9, and OcB-43, CcS-5 exhibited lower expression than CcS-16 in liver (**Figure 1D**; **Table 4A**); however, expression in none of the CcS or OcB strains exceeded the range of expression in both parental strains.

female mice of strains of strains BALB/c (*n* = 7 skin, 9 lymph nodes, 11 spleen, 9 liver), STS (6 skin, 9 lymph nodes, 8 spleen, 6 liver), CcS-5 (6 skin, 6 lymph nodes, 6 spleen, 6 liver), CcS-16 (6 skin, 6 lymph nodes, 6 spleen, 6 liver), CcS-20 (6 skin, 6 lymph nodes, 7 spleen, 6 liver), O20 (7 skin, 6 lymph nodes, 9 spleen, 6 liver), B10 (6 skin, 7 lymph nodes, 6 spleen, 6 liver), B10.O20 (5 skin, 7 lymph nodes, 6 spleen, 8 liver), OcB-9 (6 skin, 6 lymph nodes, 6 spleen, 5 liver), and OcB-43 (6 skin, 6 lymph nodes, 6 spleen, 7 liver) were compared. The data show the means ± SD. Only the differences between parental strains BALB/c and STS and strains of CcS/Dem series and parental strains O20 and B10 and strains of OcB/Dem series are shown. Nominal *P* values are shown.

Expression of *Gbp5* in skin did not differ in parental strains of CcS/Dem series BALB/c and STS (**Figure 2A**; **Table 1C**), and expression of *Gbp5* in CcS-20 exceeded the expression of both parental strains (**Figure 2A**; **Table 1C**). There was no difference in *Gbp5* expression among strains of OcB/Dem series (**Table 1C**).

We did not find any significant differences in expression of *Gbp5* among the strains of both CcS/Dem and OcB/Dem series in the lymph nodes, spleen, and liver; none of the strains in CcS/Dem or in OcB/Dem series differed from either parent (**Figures 2B–D**; **Tables 2C**, **3C**, and **4C**).

## Upregulation of *Gbp2b*/*Gbp1* and *Gbp5* mRNAs after Infection

In susceptible mice, pathology started as a nodule at the site of *L. major* infection appearing between weeks 2 and 4, which progressed in susceptible strains into a skin lesion (**Figure 3A**) (34). Strains BALB/c and CcS-16 are highly susceptible and develop large skin lesions (**Figures 3A,B**), B10.O20 develops moderate skin lesions (**Figures 3A,C**), CcS-20 is intermediate (**Figures 3A,D**) (34), and STS, CcS-5, O20, B10, OcB-9, and OcB-43 are resistant with no or minimal skin lesions (**Figure 3A**). All tested strains contain parasites in skin (**Figure 3D**) and spleen (**Figure 3E**), although the parasite load in resistant strains STS, CcS-5 and O20 (skin and spleen), and OcB-9 (spleen) is low.

Infection increased the expression of *Gbp2b*/*Gbp1* and *Gbp5* in organs of tested mice, the highest increase was observed in skin (**Figures 4**–**7**). All tested strains except CcS-5 and OcB-9 exhibited significantly higher expression of *Gbp2b*/*Gbp1* and *Gbp5* in skin after infection, irrespective of their susceptibility or resistance status (**Figure 4A**). Similarly as in uninfected mice,

TABLE 1 | Comparison of expression of *Gbp2b*/*Gbp1* or *Gbp5* among mouse strains of CcS/Dem and OcB/Dem series in skin.

*A. Differences of Gbp2b/Gbp1 expression in uninfected skin; B. Differences of Gbp2b/Gbp1 expression in skin after 8 weeks of infection; C. Differences of Gbp5 expression in uninfected skin; D. Differences of Gbp5 expression in skin after 8 weeks of infection. Green: nominal (uncorrected) P value* < *0.05; red: difference significant after correction for multiple testing at P* < *0.05.*

levels of *Gbp2b*/*Gbp1* mRNA in CcS-16 and CcS-20 exceeded those in both parental strains BALB/c and STS (Figure S1A in Supplementary Material; **Table 1B**); *Gbp5* expression in infected CcS-20 also exceeded that in both BALB/c and STS (Figure S2A in Supplementary Material; **Table 1D**).

Infection also induced increase of *Gbp2b*/*Gbp1* in inguinal lymph nodes of all strains except BALB/c and CcS-20, the highest expression was observed in CcS-5 (**Figure 5A**), which differed from all tested strains except STS (Figure S1B in Supplementary Material; **Table 2B**), but only increase of expression of B10.O20 and OcB-43 was significant after correction for multiple testing; we did not observe significant increase of *Gbp5* mRNA in lymph nodes (**Figure 5B**).

Four strains (BALB/c, STS, CcS-5, and CcS-16) show significantly increased expression of *Gbp2b*/*Gbp1* in spleen (**Figure 6A**). Expression of *Gbp5* was increased in CcS-5 (**Figure 6B**).

In liver, infection induced significant increases of *Gbp2b*/*Gbp1* mRNA in strains of CcS/Dem series, CcS-5, and CcS-16 (**Figure 7A**). Level of *Gbp2b*/*Gbp1* mRNA in CcS-16 is highest from all tested strains (Figure S1 in Supplementary Material; **Table 4B**). Expression of *Gbp5* was significantly increased in CcS/ Dem strains CcS-5 and CcS-16 and decreased in OcB/Dem strain OcB-43 (**Figure 7B**; **Table 4D**).

#### Gbp2b/Gbp1 Protein Tends to Co-Localize with *Leishmania* Parasites in Skin of Resistant and Intermediate Strains but Not in the Highly Susceptible Strain BALB/c

Expression of *Gbp2b*/*Gbp1* mRNA was highest in skin of infected mice (Figure S1 in Supplementary Material; **Figure 4**), we have therefore analyzed by immunohistochemistry a presence of Gbp2b/Gbp1 protein in the skin of selected strains BALB/c, STS, CcS-5, CcS-20, and O20 and its relationship to *L. major* parasite in infected mice. **Figure 8** shows the presence of Gbp2b/Gbp1 protein in the skin of uninfected strains. The comparison of position of *L. major* and Gbp2b/Gbp1 in the skin of chronically infected highly susceptible strain BALB/c showed few Gbp2b/Gbp1 in the vicinity of *L. major* parasites, but a large part of parasites was free of Gbp2b/

*A. Differences of Gbp2b/Gbp1 expression in uninfected lymph nodes; B. Differences of Gbp2b/Gbp1 expression in lymph nodes after 8 weeks of infection; C. Differences of Gbp5 expression in uninfected lymph nodes; D. Differences of Gbp5 expression in lymph nodes after 8 weeks of infection. Green: nominal (uncorrected) P value* < *0.05; red: difference significant after correction for multiple testing at P* < *0.05.*

Gbp1 (**Figure 9A**); the comparison of parasite load in the skin of the tested strains is shown in Figure S3 in Supplementary Material. In resistant strains STS (**Figure 9B**), CcS-5 (**Figure 9C**), and O20 (**Figure 9E**) and in intermediate strain CcS-20 (**Figure 9D**), Gbp2b/ Gbp1 co-localized with clusters of parasites (**Figures 9B–E**) that in some places formed large clusters or long stretches. Gbp2b/ Gbp1 either surrounded these clusters (**Figures 9B–D**) or formed stretches consisting of *L. major* parasites and Gbp2b/Gbp1 (**Figures 9C,E**). The tightest co-localization was observed in strains CcS-20 (**Figure 9D**) and O20 (**Figure 9E**).

#### DISCUSSION

#### Genetic Influence on Expression of *Gbp2b*/*Gbp1* and *Gbp5*

Tested strains exhibited genetic differences in *Gbps* expression both before and after *L. major* infection (**Figures 1**, **2**, **4** and **7**; **Tables 1**–**4**). Our study extends analysis of genetic influence by Staeheli and coworkers on Gbp2b/Gbp1 expression (39), who injected forty six mouse strains by poly(I);poly(C) in order to induce interferon production and tested their spleen cells for guanylate-binding activity. Tested strains were divided into Gbp2b/Gbp1 inducible and Gbp2b/Gbp1 noninducible groups. BALB/c was in the inducible group, whereas STS, O20, and C57BL/6J belonged to noninducible one (39). Our data confirm strong genetic influence on expression of *Gbp2b/Gbp1*; however, a direct comparison of outcome of study of Staeheli et al. (39) with our results is impossible due to different experimental designs. They induced Gbp2b/Gbp1 expression by poly(I);poly(C) that is structurally similar to double-stranded RNA present in some viruses, whereas we stimulated *Gbp2b*/*Gbp1* expression by the chronic infection with parasite *L. major*.

# Comparison of Genotypes in *Gbp* Cluster on Mouse Chromosome 3 Indicates *Trans*-Regulation

Our data surprisingly showed that in several organs expressions levels of *Gbp*s in recombinant congenic strains were outside the range of their parents. In skin of uninfected mice, expression

TABLE 3 | Comparison of expression of *Gbp2b*/*Gbp1* or *Gbp5* among mouse strains of CcS/Dem and OcB/Dem series in spleen.

*A. Differences of Gbp2b/Gbp1 expression in uninfected spleen; B. Differences of Gbp2b/Gbp1 expression in spleen after 8 weeks of infection; C. Differences of Gbp5 expression in uninfected spleen; D. Differences of Gbp5 expression in spleen after 8 weeks of infection. Green: nominal (uncorrected) P value* < *0.05; red: difference significant after correction for multiple testing at P* < *0.05.*

of *Gbp2b*/*Gbp1* in CcS-5, CcS-16, and CcS-20 exceeded those of both their parents BALB/c and STS (**Figure 1A**) and expression of *Gbp2b*/*Gbp1* in B10.O20 exceeded expression in parental strain B10 (**Figure 1A**). Such pattern of inheritance has been considered to be caused by *trans*-regulatory effects of non-linked or distant genes (40). The differences between parental strains and CcS/Dem strain CcS-20 persist after *L. major* infection, whereas the differences between expression of parents and CcS-5 and CcS-16 and between parent B10 and the strain B10.O20 disappear after infection (**Figure 1A**; Figure S1A in Supplementary Material; **Tables 1A,B**). Expression of *Gbp5* in skin of uninfected CcS-20 exceeded level of both parents (**Figure 2A**; **Table 1C**) but was significantly higher only than the parental strain STS after 8 weeks of infection (Figure S2A in Supplementary Material; **Table 1D**). CcS-5 and CcS-16 highly differed in the expression of both *Gbp1*/*Gbp2b* and *Gbp5* in lymph nodes and liver of infected mice; these strains also differed in expression of *Gbp5* in spleen (**Tables 2B,D**, **3D**, and **Tables 4B,D**).

Comparison of genotypes of the tested strains (33, 41, 42) (this study) in the *Gbp* cluster on the mouse chromosome 3 (**Figure 10**) revealed that strains CcS-5, CcS-16, and CcS-20 exhibiting higher expression of *Gbp* had *Gbp* genotype identical to that of BALB/c (C). Similarly, highly differing CcS-5 and CcS-16 strains carry the same *Gbp* allele. The presence of the same allele of *Gbp* gene cluster on chromosome 3 in strains that differ in other genes suggests that their differences in expression of *Gbp2*/*Gbp1* and/or *Gbp5* from one or both parents or from other RC strain are due to regulatory influence of non-*Gbp* gene(s) of STS origin carried on other genetic segments (*trans*-regulation). In the OcB/Dem series, B10.O20 carried in *Gbp* cluster B10 genotype (B), which similarly indicated *trans*-regulation of expression from O20 genome situated outside *Gbp* cluster (**Table 1A**; **Figure 10**). This *trans*-regulation can be partly overlaid by other regulatory events appearing after infection. Further genetic studies will be needed to elucidate nature of regulatory events observed in our studies.

The observations of progeny having a phenotype, which is beyond the range of the phenotype of its parents, are not rare. For example, analysis of gene expression from livers in chromosome substitution mouse strains revealed that only 438 of the 4,209 expressed genes were inside the parental range (40).

TABLE 4 | Comparison of expression of *Gbp2b*/*Gbp1* or *Gbp5* among mouse of CcS/Dem and OcB/Dem series strains in liver.

*A. Differences of Gbp2b/Gbp1 expression in uninfected liver; B. Differences of Gbp2b/Gbp1 expression in liver after 8 weeks of infection; C. Differences of Gbp5 expression in uninfected liver; D. Differences of Gbp5 expression in liver after 8 weeks of infection. Green: nominal (uncorrected) P value* < *0.05; red: difference significant after correction for multiple testing at P* < *0.05.*

These observations are due to multiple regulatory interactions, which in new combinations of these genes in recombinant congenic or chromosomal substitution strains can lead to the appearance of new phenotypes that exceed their range in parental strains.

#### Increased Expression of *Gbp2b*/*Gbp1* and/or *Gbp5* in Resistant Mice Suggests Hidden Inflammation

We and others have demonstrated that *Leishmania* parasites are present not only in organs of infected susceptible mice with clinical manifestations of the disease but also in clinically asymptomatic mice of resistant strains (37, 43–46). This is also shown in **Figures 3D,E** and **9B,C,E**. **Figures 4**–**7** show that the expression of *Gbp2b*/*Gbp1* and/or *Gbp5* has increased after infection in at least one organ of each of the tested mice, including the resistant ones (STS, CcS-5, O20, B10, OcB-9, and OcB-43), which had no or only minimal and transient clinical symptoms. This strongly suggests that persistent parasites can contribute to the maintenance of protective immunity, which was manifested in our experiments by the increased levels of *Gbp2b/Gbp1* and *Gbp5* in resistant mice. It was demonstrated previously that this latent infection is controlled by inducible nitric oxide synthase (43) and phagocyte NADPH oxidase (46). It remains to be established, whether defense mechanisms including Gbps that were found to act against other pathogens (16, 23), operate also in *Leishmania*infected mammalian host. In defense against *M. bovis*, Gbp2b/ Gbp1 and Gbp7 could promote NADPH oxidase activity after the recruitment of gp91phox and gp22phox components to bacteria vacuoles (23), whereas parasite *T. gondii* was directly attacked via Gbp supramolecular complexes (16). The observed association (**Figure 9**) of Gbp2b/Gbp1 with *L. major* parasites in the skin of resistant and intermediate strains but not the highly susceptible strain BALB/c may suggest a role of this protein in response against the *L. major* pathogens.

Importantly, persistent parasites, besides stimulating protective immune reactions, can also represent a danger for hosts (45). The increased expression of *Gbp2b*/*Gbp1* and *Gbp5* in clinically asymptomatic mice reveals the price exacted from the

highly susceptible strains) till week 8 (end of experiment) (A). Median values of skin lesions of mice tested in expression experiments are shown. Skin lesion caused by *Leishmania major* in female mouse of BALB/c strain at week 8 after infection (B). Skin lesion caused by *L. major* in female mouse of B10.O20 strain at week 8 after infection (C). Parasite load in skin (D) and spleen (E) of infected female mice of strains BALB/c (*n* = 5 skin, 11 spleen), STS (13 skin, 7 spleen), CcS-5 (11 skin, 4 spleen), CcS-16 (6 skin, 6 spleen), CcS-20 (7 skin, 9 spleen), O20 (6 skin, 12 spleen), B10 (9 skin, 13 spleen), B10.O20 (7 skin, 15 spleen), and OcB-9 (7 skin, 7 spleen). The data show the means ± SD.

FIGURE 5 | Differences in expression of *Gbp2b*/*Gbp1* and *Gbp5* in inguinal lymph nodes of uninfected and infected mice. Expression of *Gbp2b/Gbp1* (A) and *Gbp5* (B) in inguinal lymph nodes of uninfected and infected female mice of strains BALB/c (11 infected and 9 uninfected), STS (6 infected and 6 uninfected), CcS-5 (6 infected and 6 uninfected), CcS-16 (6 infected and 6 uninfected), CcS-20 (6 infected and 6 uninfected), O20 (9 infected and 6 uninfected), B10 (13 infected and 7 uninfected), B10.O20 (11 infected and 7 uninfected), OcB-9 (7 infected and 6 uninfected), and OcB-43 (5 infected and 6 uninfected) were compared. Animals were subcutaneously inoculated with 107 promastigotes of *Leishmania major.* Control, uninfected mice were kept in the same animal facility. Both groups were killed after 8 weeks of infection. The data show the means ± SD from 12 independent experiments. Nominal *P* values are shown.

of uninfected and infected female mice of strains BALB/c (14 infected and 11 uninfected), STS (12 infected and 8 uninfected), CcS-5 (6 infected and 6 uninfected), CcS-16 (6 infected and 6 uninfected), CcS-20 (5 infected and 7 uninfected), O20 (9 infected and 9 uninfected), B10 (10 infected and 6 uninfected), B10.O20 (6 infected and 6 uninfected), OcB-9 (6 infected and 6 uninfected), and OcB-43 (6 infected and 6 uninfected) were compared. Animals were subcutaneously inoculated with 107 promastigotes of *Leishmania major.* Control, uninfected mice were kept in the same animal facility. Both groups were killed after 8 weeks of infection. The data show the means ± SD from 12 independent experiments. Nominal *P* values are shown.

organism by a dormant infection. This finding deserves attention, because elevated levels of human GBP1 are directly involved in the endothelial dysfunction and the regulation of endothelial progenitor cells activity in patients with the autoimmune diseases such as rheumatoid arthritis, systemic lupus erythematosus and systemic sclerosis (47). In mice, elevated levels of Gbp3 and Gbp6 were linked with the pathogenesis of atherosclerosis (48). In humans with colorectal cancer, the anti-angiogenic effect of increased levels of GPB1 was beneficial in colorectal carcinoma patients, where it was associated with sustained reduction of intratumoral angiogenic activity and improved cancer-related survival (49).

FIGURE 7 | Differences in expression of *Gbp2b*/*Gbp1* and *Gbp5* in liver of uninfected and infected mice. Expression of *Gbp2b/Gbp1* (A) and *Gbp5* (B) in liver uninfected and infected female mice of strains BALB/c (13 infected and 9 uninfected), STS (6 infected and 6 uninfected), CcS-5 (6 infected and 6 uninfected), CcS-16 (6 infected and 6 uninfected), CcS-20 (7 infected and 6 uninfected), O20 (8 infected and 6 uninfected), B10 (12 infected and 6 uninfected), B10.O20 (11 infected and 8 uninfected), OcB-9 (7 infected and 5 uninfected), and OcB-43 (7 infected and 6 uninfected) were compared. Animals were subcutaneously inoculated with 107 promastigotes of *Leishmania major.* Control, uninfected mice were kept in the same animal facility. Both groups were killed after 8 weeks of infection. The data show the means ± SD from 12 independent experiments. Nominal *P* values are shown.

(C), CcS-20 (D), and O20 (E) mice infected for 8 weeks with *L. major* were stained with the anti-*Leishmania* lipophosphoglycan mouse monoclonal antibody (cat. no. CLP003A, Cedarlane, Hornby, Canada) and TRITC labeled IgM (115-025-020, Jackson ImmunoResearch, West Grove, PA) all diluted 1:500 and the rabbit anti-Gbp1 Polyclonal antibody (PA5-23509, Thermo Fisher Scientific, Rockford, IL, USA) diluted 1:100 and anti-rabbit-AlexaFluor-647 (cat. no. 711-605-152; Jackson ImmunoResearch, West Grove, PA) diluted 1:500. Nuclei of the cells were stained with bisBenzimide H33258 (Sigma-Aldrich, St. Louis, MO, USA) 10 mg per 1 ml diluted 1:1,000. Images were captured with microscope Leica DM6000 objective HCX PL Apo 40×/0.75 PH2 and color camera Leica DFC490. Evaluation of images was done with Fiji ImageJ 1.51n. Figures are representatives of data from 8 to 11 mice (see Materials and Methods) in 3 of them 10 fields (320.66 × 239.57 µM) from each mouse were analyzed, in the rest one field was analyzed to verify the results. White arrow shows *L. major* amastigotes (red color), green arrows show Gbp2b/Gbp1 protein (yellow color), red arrows point to amastigotes co-localized with Gbp2b/Gbp1, whereas yellow arrows show either Gbp2b/Gbp1 surrounding parasite clusters or stretch of parasites and Gbp2b/Gbp1. Cell nuclei are stained in blue.

The immune reactions accompanying persistent *Leishmania* infection might be very important, because in addition to 12 million people presently suffering from the clinical manifestations of leishmaniasis (50), there are at least 120 million people with asymptomatic infection (45). It needs to be elucidated, whether such clinically asymptomatic people harboring persistent *Leishmania* parasites are more prone to immune-related diseases.

## CONCLUSION

Our results represent the presently most comprehensive information about expression of *Gbps* in leishmaniasis *in vivo*.

We found that expression of *Gbp2b*/*Gbp1* and *Gbp5* is under strong genetic control involving in some strains also *trans*regulation both in uninfected and *L. major*-infected mice.

We have observed that in several organs, expression of *Gbps* in recombinant congenic strains was outside the range of their parents. Tests of different strains that carry the same *Gbp* cluster genotypes on chromosome 3 indicate a *trans*-regulation of *Gbp2b*/*Gbp1* and *Gbp5* by genes that are not closely linked to *Gbp* genes. This finding may open way to identification and manipulation of these presently unknown genes.

Our results also point out that expression of *Gbp2b*/*Gbp1* and *Gbp5* was increased even in organs of resistant mice, which might suggest a hidden inflammation. It remains to be established whether the clinically asymptomatic infection might represent danger in predisposing organism to other diseases.

Co-localization of Gbp2b/Gbp1 protein with most *L. major* parasites in skin of resistant and intermediate strains STS, CcS-5, O20, and CcS-20 but not in highly susceptible BALB/c mice suggests that this molecule might play role in defense against leishmaniasis and opens new research direction in analysis of control of persistent parasites.

#### ETHICS STATEMENT

All experimental protocols utilized in this study comply with the Czech Government Requirements under the Policy of

#### REFERENCES

Animal Protection Law (No. 246/1992) and with the regulations of the Ministry of Agriculture of the Czech Republic (No. 207/2004), which are in agreement with all relevant European Union guidelines for work with animals and were approved by the Institutional Animal Care Committee of the Institute of Molecular Genetics AS CR and by Departmental Expert Committee for the Approval of Projects of Experiments on Animals of the Academy of Sciences of the Czech Republic (permissions Nr. 190/2010; 232/2012).

#### AUTHOR CONTRIBUTIONS

YS planned and performed parasitology and expression experiments and contributed to the writing of the manuscript. VV performed parasitology experiments, analyzed the data, and contributed to the writing of the manuscript. TK contributed to the writing of the manuscript, parasitology experiments, estimation of parasite numbers, and data analysis. HH designed and performed immunohistochemistry analysis and analyzed the data. IK contributed to the estimation of parasite numbers. MS performed parasitology experiments. PD analyzed the data and contributed to the writing of the manuscript. ML conceived the study, interpreted data, and wrote the manuscript. All authors reviewed the manuscript.

#### ACKNOWLEDGMENTS

YS and IK are PhD students of Faculty of Science, Charles University, Prague, Czech Republic. MS is a PhD student of Faculty of Science, The University of South Bohemia, České Budějovice, Czech Republic. We thank J. Bartůněk and T. Pokorná for technical help.

#### FUNDING

This work was funded by the Czech Science Foundation (https:// gacr.cz/en/) (Grants GACR 14-30186S, GACR 16-22346S, and 13-41002P). It was also supported by the Czech Academy of Sciences (http://www.avcr.cz/en/) (Project grant RVO: 68378050).

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at http://www.frontiersin.org/articles/10.3389/fimmu.2018.00130/ full#supplementary-material.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2018 Sohrabi, Volkova, Kobets, Havelková, Krayem, Slapničková, Demant and Lipoldová. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*


**5. Antifungal assay**

The maximum extent and rate of activity is known to vary with the seed culture mediums, the physiological age of the culture, and the type of culture medium. The initial assay was achieved using two-fold serial dilution method against *Saccharomyces cerevisiae* ATCC 7754. This yeast was purchased from American Type Culture Collection (Manassas, VA).


*Article* **Comparative Study of Methods for Cycle Length Estimation in Fractionated Electrograms of Atrial Fibrillation**

**Diego Osorio 1, Aikaterini Vraka 1, José Moreno-Arribas 2, Vicente Bertomeu-González 2, Raúl Alcaraz 3 and José J. Rieta 1,\***

**Abstract:** Atrial cycle length (CL) is an important feature for the analysis of electrogram (EGM) characteristics acquired during catheter ablation (CA) of atrial fibrillation (AF), the commonest cardiac arrhythmia. Nevertheless, a robust ACL estimator requires the precise detection of local activation waves (LAWs), which still remains a challenge. This work aims to compare the performance in (CL) estimation, especially under fractionated EGMs, of three different LAW detection methods relying on different operation strategies. The methods are based on the hyperbolic tangent (HT) function, an adaptive amplitude threshold (AAT) and a (CL) iteration (ACLI), respectively. For each method, LAW detection has been assessed with respect to manual annotations made by two experts and performance has been estimated by confusion matrix and mean and individual (CL) error calculation by EGM types of fractionation. The influence of EGM length on the individual (CL) error has been additionally considered. For the HT method, accuracy, sensitivity and precision were 92.77–100%, while for the AAT and ACLI methods they were 78.89–99.91% for all EGM types. The CL error on the HT method was lower than AAT and ACLI methods (up to 12 ms versus up to 20 ms), with the difference being more prominent in complex EGMs. The HT method also showed the lowest dependency on EGM length, presenting the lowest and least variable error values. Therefore, the HT method achieves higher performance in (CL) estimation in comparison with previous LAW detection techniques. The high robustness and precision demonstrated by this method suggest its implementation on CA mapping devices for a more successful location of ablation targets and improved results during CA procedures.

**Keywords:** atrial fibrillation; electrogram; complex fractionated atrial electrograms; local activation waves; detection; invasive recordings; cycle length; comparison

#### **1. Introduction**

Being the most common cardiac arrhythmia in clinical practice, atrial fibrillation (AF) affects more than 43 million people worldwide [1]. The risks of stroke, myocardial infarction, heart failure or dementia are increased in patients suffering from AF, which is also associated with increased mortality [1,2]. The ever-growing life expectancy in comparison to the high AF prevalence for older individuals [1] and the significant economic burden that AF provokes [1,3] highlight the need for an efficient treatment.

So far, catheter ablation (CA) is considered the star AF treatment due to its high effectiveness and safety [1,4,5]. The aim of the CA procedure is to eliminate the ectopic atrial beats that can trigger or sustain AF by electrically isolating the areas that are thought to present arrhythmogenic activity or fibrosis, respectively [1,6]. Pulmonary veins (PVs) are the most common AF foci, especially for patients in paroxysmal AF [1,7,8]. Nevertheless, other atrial sites may contribute to AF perpetuation, especially in persistent AF patients,

**Citation:** Osorio, D.; Vraka, A.; Moreno-Arribas, J.; Bertomeu-González, V.; Alcaraz, R.; Rieta, J.J. Comparative Study of Methods for Cycle Length Estimation in Fractionated Electrograms of Atrial Fibrillation. *J. Pers. Med.* **2022**, *12*, 1712. https://doi.org/10.3390/ jpm12101712

Academic Editor: José Miguel Rivera-Caravaca

Received: 25 August 2022 Accepted: 11 October 2022 Published: 13 October 2022

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

where the atrial remodeling takes place to a higher extent [1,9]. Consequently, additional CA may be necessary in order to achieve freedom from AF [10]. Areas targeted for additional CA are defined via electrogram (EGM) analysis, which allows the detection of remodeled tissue favoring AF perpetuation or even triggering AF. Low-voltage, high dominant frequency (DF), low AF cycle length (CL) or highly complex EGMs, called complex fractionated atrial electrograms (CFAEs), are some of the most established techniques to define candidate CA targets, albeit showing controversial results [11–16]. These types of EGMs are thought to correspond to areas of slow conduction or pivotal points of reentrant wavelets, revealing the presence of atrial fibrosis which provokes AF perpetuation [17].

For many of the aforementioned techniques, the detection of local activation waves (LAWs) is an elemental step. In low fractionated EGMs, LAWs can be easily detected, hence facilitating the analysis [18]. However, in CFAEs and highly fractionated EGMs, LAWs detection can become a hard task often misleading the CL estimation and, consequently, the procedures and assessments based on this feature [13]. Moreover, the possibility of ambiguous LAWs annotations in CFAEs may be an explanation for the low impact of non-PV CA on the freedom from AF [14,15]. Therefore, the development and disposal of a reliable and robust LAW estimator are of paramount importance. Although various approaches to the development of LAWs detectors have been attempted so far [19–27], most of them focus on unipolar EGMs [20,21,25–27], where LAWs detection is a much easier task [13]. Nevertheless, as LAWs detection in bipolar EGMs is more complicated due to the nonstationary nature of the recorded dynamics as well as the high dependency on the wavefront direction, not but a few techniques for bipolar EGMs have been developed so far [19,22–24].

One of the first and most cited works employed an adaptive amplitude-based threshold which varied according to the last detected peaks, showing low error values for regular or low fractionated EGMs but increased rates for fractionated signals [19]. The CL-based method is another highlighted work [22]. Search is performed with an initial amplitude threshold which is then modified according to CL. This method has also shown low error rates for regular EGMs. Nevertheless, when CL becomes irregular, as normally happens in CFAEs, error rates increase. The DF-based method detects LAWs corresponding to periodic peaks of the EGMs [23]. Although this method showed good performance compared to previous methods regarding periodic EGMs, it also failed to show satisfactory results for aperiodic signals, which is a much more frequent event in highly fractionated EGMs. Thus, it will not be considered in this study. Finally, the newer of the aforementioned methods is based on the application of the hyperbolic tangent (HT) function on CFAEs in order to highlight and facilitate low amplitude detections as well as to limit very high amplitude activations using a hybrid methodology combining attributes of the CL- and adaptive amplitude-methods, showing high performance in CFAEs and complex EGMs [24].

Given the advantage of bipolar over unipolar EGMs in being insensitive to ventricular activity and the high use of the former in AF studies, the appropriate strategy choice for works that may need to implement LAWs detection in their analysis can be quite confusing yet inevitable. Despite the results shown in the aforementioned studies, each one of them has used a different database. Hence, a direct comparison would be unfair or would not be directly performed. The main objective of the present work is to perform a straightforward comparison between the three most robust aforementioned LAW detectors, the adaptive amplitude threshold (AAT) [19], the atrial CL iterator (ACLI) method [22] and the hyperbolic tangent (HT) method [24], by reproducing each one of them and comparing the results over the same database. Since the common denominator of LAWs detection in any of the aforementioned studies is the CL estimation, using this parameter as a performance evaluation metric seems a fair option. Comparison is performed over a range of factors that could affect or create biased results, such as the length of EGMs or different fractionation levels. In this way, a multidimensional assessment can be performed in order to offer a more reliable and direct comparative measure for future studies.

The manuscript is organized as follows. Section 2 presents the dataset used in this study and the preprocessing methods applied, as well as a brief presentation of the methodology of the strategies under comparison and the evaluation methods used. Section 3 includes the results of the comparison, while Section 4 provides an overview of the methods under comparison and analyzes the main conclusions derived from the results, which are then summarized in Section 5.

#### **2. Methods**

The dataset employed in this study consisted of 119 EGMs of 10 s in length, recorded by a Labsystem™ PRO EP recording system (Boston Scientific, Marlborough, MA, USA) after obtaining the written consent of 22 patients undergoing CA procedures. EGMs were sampled at 1 kHz and filtered by an adaptive notch filter to reduce powerline interference together with a 0.5–500 Hz band-pass filter to limit high-frequency noise. The dataset was manually annotated and classified by EGM types by two expert physicians according to Wells' classification [28]. Type I EGMs consisted of clear LAWs and almost isoelectric baseline, type II contained low perturbations in the baseline and a higher number of deflections at each LAW while type III contained a high fluctuating baseline and lowamplitude fragmented LAWs. Sixteen EGMs were classified by the experts as AF type I, 19 as AF type II and 84 as AF type III, so that the vast majority of the EGMs selected pertained to the most difficult AF type to delineate. Signals were denoised by a Wavelet Transform-based technique, which has been proved to perform better than regular filtering, providing quick results, effective noise elimination and optimal preservation of the signal morphology [29].

#### *2.1. LAW Detection Methods*

The methodology of each of the studied LAW detection strategies is presented briefly in the next subsections, followed by a comparative example of the preprocessing method adopted by each strategy.

#### 2.1.1. Adaptive Amplitude Threshold Method

The first method implemented was the AAT algorithm [19]. This technique is based on Botteron and Smith preprocessing [18], a technique developed to generate waveforms proportional to the amplitude of the EGM components with frequencies within the bandpass filter cutoffs. This method includes 40–250 Hz band-pass filtering, followed by signal rectification and low-pass filtering with the cut-off frequency at 20 Hz. The AAT algorithm employs an AAT based on the last ten detected activations with exponentially decreasing weights (*sw*), resulting from the application of the Botteron filtering, using a blanking period of 55 ms between activations. Lastly, the modulus of the original signal is filtered by a non-causal moving average filter with 90 coefficients, resulting in a new signal *sf* , where the positive zero crossings closer to a local peak of *sw* define the activation times. A more detailed description of the algorithm can be found elsewhere [19].

#### 2.1.2. Atrial Cycle Length Iteration Method

The ACLI method [22] is based on a modified Botteron and Smith [18] preprocessing, replacing the band-pass filter with a 40 Hz high-pass filter and increasing the low-pass filter cut-off to 30 Hz. The algorithm performs an iterative LAW detection, starting with the highest amplitude peak and moving to the next peak according to an amplitude descending order, making use of a blanking period of 50 ms, until the average CL is lower than 275 ms. Once this condition is fulfilled, a second condition has to be accomplished so that the iteration stops: either the mean CL has to be less than the median CL plus 5 ms or the amplitude of the current peak has to be 20% less than the amplitude of the previous peak. Finally, a loop control checks the existence of intervals longer than 1.5 × the median CL and adds, in case such intervals are found, the highest peaks within these intervals to the activation set. More details on the methodology can be found elsewhere [22].

#### 2.1.3. Hyperbolic Tangent Method

The HT method [24] also modifies the Botteron and Smith preprocessing steps, expanding the range of the bandpass filter cutoffs from 40–250 to 20–250 Hz. Studies on the original preprocessing effect on CFAEs demonstrated that slower components of activations were being lost, provoking an insufficient waveform response, leading the algorithm to miss the respective activations [24]. Therefore, a reduction in the minimum band-pass filter cut-off was necessary, with 20 Hz being the optimal frequency based on analyses performed.

The preprocessing of this method also used a novel technique, the HT function, aiming to reduce high amplitude variability observed in CFAEs in order to facilitate the detection process that follows next. For this purpose, each recording is segmented into 7-second intervals with 25% overlapping and classified according to their fractionation level by their kurtosis calculated at each 1 s and averaged, so that the HT function is applied exclusively on CFAEs. The detection process is initiated by an amplitude-based localization of all peaks higher than 0.4 mV, taking into consideration the minimum refractory period of 50 ms. The process continues with the median CL calculation and a secondary search on intervals longer than the median CL is performed, reducing the amplitude threshold proportionally to the difference between the median CL and the current interval. The process finishes with the calculation of the barycenter of each LAW, defined as the sample point that equally divides the area of the modulus of the LAWs, in order to estimate the activation time. More details on this algorithm are provided elsewhere [24].

A comparative example of these three different preprocessing strategies applied over the same EGM is shown in Figure 1, where it can be observed how the HT method produces a signal in which activations are easier to detect automatically in later stages because small peaks are enhanced and large peaks are limited. Additionally, block diagrams of each technique's methodology are presented in Figure 2.

**Figure 1.** Comparative analysis of the preprocessing stage of each of the compared LAW detection methods applied on the same EGM: (**a**) Denoised EGM to be processed; (**b**) Result of AAT preprocessing; (**c**) Result of ACLI preprocessing; (**d**) Result of HT preprocessing. The dashed line in (**d**) corresponds to the application of the hyperbolic tangent.

**Figure 2.** Block diagrams of the LAW detection methods compared. Underlined values in the first blocks show modified thresholds with respect to traditional Botteron's preprocessing method. BP: band-pass filtering; LP: low-pass filtering; HP: high-pass filtering; *Pi*: amplitude of the *i*-th peak; med{CL}: median cycle length.

#### *2.2. Performance Evaluation*

Confusion matrix and CL error metrics were employed in order to evaluate the performance. Regarding the confusion matrix, *Accuracy*, *Sensitivity* and *Precision* were calculated. True positive (*TP*) were the LAWs annotated both by the experts and by the method under analysis. False positives (*FP*) were the LAWs that were only annotated by the method, while false negatives (*FN*) were the LAWs that were annotated only by the experts. Since there are no true negatives (*TN*), these metrics were calculated as follows:

$$Accuracy = \frac{TP}{TP + FP + FN} \,\prime \tag{1}$$

$$Sensitivity = \frac{TP}{TP + FN} \tag{2}$$

and

$$Precision = \frac{TP}{TP + FP}.\tag{3}$$

Details on CL error calculation are provided next. The effect of EGMs duration as well as fractionation has also been considered.

#### 2.2.1. Detection Performance

The number of LAWs successfully detected at each method was calculated by comparing the annotations made by the algorithm with the manual annotations made by the experts, allowing a threshold of 40 ms. Performance was evaluated in terms of Accuracy, Sensitivity and Precision, providing an overall perspective on performance, including the undersensing and oversensing rates of each algorithm.

#### 2.2.2. Mean CL Error Estimation

The error in estimating the average CL at any of the three methods was measured as the absolute difference between the mean CL of the manual and algorithm's annotations, as can be seen from Figure 3, which will be explained more in-depth later for individual CL error estimation. Nevertheless, considering that oversensing or undersensing can have an effect on the mean CL, producing misleading results as can be seen in Figure 3B, mean CL error can not be considered as a reliable performance indicator. Thus, an analysis of individual CL errors was additionally developed in order to obtain more realistic results.

**Figure 3.** Comparative example for calculating the CL estimation error. (**A**) EGM with experts' (circle) and method's annotations: (a) expert's CL (b) method's CL and (c) errors in CL estimation. Top: Situation 1 (squares) undersensing case (missing activation-red). The arrow shows the starting points in case of unequal CLs for the undersensing case; Bottom: Situation 2 (triangles) undersensing and oversensing case (false positive activation-red). (**B**) Error estimation for the two example situations. (left) Mean CL error; (right) Individual CL error. Notice how the simultaneous existence of underand oversensing in Situation 2 may create erroneous CL calculation for the mean CL case.

#### 2.2.3. Individual CL Error Estimation

As previously stated, this metric was introduced in order to reliably calculate the CL error for each algorithm, without any bias inserted from possible under- or oversensing tendencies. In this case, the error at each CL is defined by the absolute difference between the manual annotations, called reference annotations, and the algorithm's annotations. The individual CL error is finally calculated by the sum of errors found at each CL divided by the number of reference intervals. The formula for the calculation of the individual CL error is given in the right column of Figure 3B. If the CL was the same for the manual and algorithm annotations, the respective error was set to zero. However, if there was some under- or oversensing detected, the evaluation process was able to quantify the error.

For the undersensing case, as the given CL was longer for the algorithm than for the manual annotations case, the corresponding error is calculated by subtracting the reference CL from the algorithm's CL. However, in that case, the next interval under comparison will involve the following reference segment, while the algorithm's segment will remain the same as before, considered to be starting from the reference segment's initiation point so that a double counting on error will be avoided. An example of this process can be seen in Figure 3A, with experts' (circle) and method's (square, triangle) annotations. The top panel shows Situation 1 (square) with an undersensing case where the arrow shows the beginning of the algorithm's segment (missed activation-red). For each interval between annotations, line (a) indicates the expert's CL, line (b) the method's CL and line (c) errors in the estimation of each CL. The bottom panel shows Situation 2 (triangles) with both undersensing (missed activation-red) and oversensing cases (falsely detected activationred). Regarding the oversensing case, the error at the given interval will be calculated by subtracting the algorithm's CL from the reference CL, for each one of the algorithm's intervals between a reference interval, as can be seen from Situation 2 of Figure 3A. Finally, as the mean CL of Situation 1, with one missed activation, is 120 ms, then the mean error is 20 ms. Nevertheless, the simultaneous existence of under- and oversensing in Situation 2 will create an erroneous mean CL calculation of 0 ms, as this method had a mean CL of 100 ms, even though not all the activations corresponded to the experts' annotations. On the contrary, the individual error of Situation 2 was 33.33 ms, both indicating a nonperfect annotator and showing a higher error than Situation 1 (16.67 ms), which only had one mishit.

#### 2.2.4. Influence of Electrograms Duration

Another aspect analyzed in this study was the influence of the EGMs duration on the individual CL error of each method. For each type of EGM, the three methods were applied progressively from 1 to 10 s of duration over the dataset, computing the individual CL error as a measure of the CL estimation precision.

#### **3. Results**

#### *3.1. Detection Results*

An example of the three AF-type EGMs along with LAWs annotations from experts and the compared techniques can be seen in Figure 4. While LAWs detection is pretty straightforward for AF types I and II, LAWs detection in AF type III EGMs is significantly complicated, with each algorithm presenting some false positive or negative detections. Results on detection performance for the methods under comparison are shown in Table 1 for EGMs of all AF types. Both ACLI and HT methods show optimal results for types I and II. On the contrary, the AAT method showed lower accuracy and precision rates, which is especially unusual for AF type I. This fact is provoked by the high variability in amplitude observed in several type I EGMs of the dataset used for the present analysis. Figure 3 corresponds to a type I EGM example, in which the impact of amplitude variability on the AAT method's detection performance can be seen. Although AAT's accuracy rate is not lower than 90% in type I EGMs, this relatively low performance implies significant deficiencies of the AAT method and a higher CL error, provoked by the longer CL of AF type I EGMs, leading to the lower number of intervals and, consequently, the higher weight of each mishit.

**Table 1.** Electrogram detection performance results for all AF types of the three LAW detectors compared in this study in terms of Accuracy (*Acc*), Sensitivity (*Se*) and Precision (*Pr*).

The results for type III EGMs, also shown in Table 1, are of great importance in understanding the performance of each method, as this type of EGM coincides with the most fragmented signals and CFAEs. Firstly, the AAT method shows the lowest accuracy results, while sensitivity and precision are quite compensated, with values of 85.38% and 91.11%, respectively, due to higher undersensing than oversensing incidents. The ACLI method shows improved results in comparison with the AAT method. However, a quite intense discrepancy can be seen between undersensing and oversensing tendencies of this method, as can be observed from the difference between sensitivity and precision rates at 87.51% and 96.89%, respectively. Finally, the HT method presents the most satisfactory results for all confusion matrix parameters, with a balance observed between undersensing and oversensing tendencies, expressed by the 95.30% and 97.24% of sensitivity and precision values, respectively, while showing the highest accuracy rates as well.

**Figure 4.** One-second example of the three EGM types and the annotations made from the three algorithms as well as from the experts. Experts' annotations are considered as ground truth. For AF type I, activations were clearly distinguished from the baseline and all the methods were able to identify them. For AF type II, although the baseline shows little perturbation, lower amplitude activations slightly complicate the procedure, leading to a mishit for the AAT algorithm. Finally, for AF type III EGMs, mishits and false positives are observed for all methods across the EGM.

#### *3.2. Cycle Length Estimation*

Regarding CL results, Figure 5a shows the mean and individual CL errors on type I EGMs. The mean CL error in ACLI and HT methods is close to 0, as expected from the high detection accuracy and the equilibrium between sensitivity and precision. Individual CL error in these methods is also quite low, yet slightly higher than mean CL error since bias from undersensing or oversensing behavior is removed in this parameter. Regarding AAT results, both mean and individual CL error are significantly higher, at about 10 ms and 12 ms, respectively. This is explained by the relatively low accuracy of this method in type I EGMs. As before, higher individual CL errors are explained by the balance between the undersensing and oversensing rates of these EGM types.

**Figure 5.** Comparative analysis of errors estimating CL using the mean CL and the individual CL methods for the three LAW detectors compared: (**a**) Results for type I electrograms; (**b**) type II electrograms; and (**c**) type III electrograms. Values indicated with + stand for outliers.

The results for type II EGMs are presented in Figure 5b. ACLI and HT methods show results similar to those from type I EGMs, with an increased difference between individual and mean CL errors. As previously mentioned, mean CL estimation for the AAT method compensated for the errors provoked by under- or oversensing. Lower amplitude variability in the employed dataset led to better performance in AF type II than type I EGMs for the AAT method, as was observed from the detection results. According to this performance improvement, the error in type II EGMs was lower than type I for the AAT method.

Finally, Figure 5c shows the CL error on type III EGMs, which include CFAEs. Being the most complicated case, analysis of this type of EGM clearly demonstrates that individual CL error estimation is a more suitable parameter, efficiently removing any impact of external factors on the CL error calculation, as can be seen from the discrepancies between the mean and individual CL error estimation of each method, leading to invalid conclusions in the case of the AAT and ACLI methods. Error on all methods is higher in this EGM type, with the HT technique showing the best results both in mean (about 4 ms in HT vs. 9 ms in AAT vs. 11 ms in ACLI) and individual (about 12 ms in HT vs. 20 ms in AAT vs. 17 ms in ACLI) CL error analysis.

#### *3.3. EGM Length Influence*

The results of the application of each method over the dataset for different EGM lengths are presented in Figure 6, where individual CL error is analyzed. The HT method shows the most stable results regarding length, being additionally improved in comparison with the other methods. For type I EGMs, the ACLI method shows similar results. As the fractionation level increases, however, ACLI performance falls, showing significantly worse results. Finally, AAT error is the highest regarding all AF types, while this method also shows a higher dependency on the EGM length.

**Figure 6.** Individual CL error depending on EGM length for the three LAW detection methods analyzed: (**a**) Results for type I electrograms; (**b**) type II electrograms; and (**c**) type III electrograms.

#### **4. Discussion**

Due to the high importance of precise AF substrate mapping, several methods for LAW detection have been developed [19,20,22–26] until the present. However, some of these techniques cannot be reproduced or directly applied on CA devices, due to either the particularity of the database [20] used or to the controversial methodology followed, using computational models premised on assumptions about tissue's geometry [25,26]. Other techniques can be directly implemented on substrate mapping devices, such as the AAT [19], the ACLI [22], the dominant frequency (DF) [23] and the HT [24] methods. Each one of these studies has shown promising results, the correct choice for future works in need of a reliable LAWs detector, though, can be confusing as the database used is different.

The present study has performed a direct comparison of three of the aforementioned methods, the AAT, the ACLI and the HT method. The reason why the DF method was not included in the comparison is the lack of individual CL estimation of the DF method, which is one of the main evaluation parameters of the present study. In addition, the DF method failed to provide satisfactory results under aperiodic signals, which is the most typical situation in highly fractionated EGMs.The results indicate that HT is the most robust and reliable method, showing higher accuracy and precision and lower CL errors, suggesting its recruitment for both future studies and CA devices.

The high performance of the HT method can be assigned to the alternative preprocessing strategy and the recruitment of the HT function, which prevented the detector from missing low amplitude LAWs in CFAEs. What is more, this method also adopted a segmentation strategy that allowed the adaptation to various signal lengths, while the combination of an adaptive amplitude threshold with CL correction allowed the detection of LAWs in segments with high amplitude variability [24]. Finally, activation time estimation through the barycenters led to an estimation closer to the experts' annotation.

#### *4.1. Evaluation on Different Levels of Fractionation*

EGM fractionation plays a significant role in the performance of each method. Normally, in low fractionation environments, performance is high and only a few errors are observed. Thus, low performance in low fractionation EGMs implies a lack of reliability. As EGMs become more fractionated, with indiscernible activations, the LAW detector's accuracy tends to drop. However, since CFAEs, which are highly fractionated EGMs, are the main challenge of any LAW detector, the performance of each method on that kind of recording is crucial and serves as an index to evaluate reliability and robustness.

This study has utilized the Wells' classification for types of EGMs, with ascending levels of fractionation from type I to type III EGMs. In this way, insufficient results in type I EGMs indicate a less reliable method, while competent results on type III EGMs affirm the reliability of another method. The final assessment is a trade-off among the performance of each EGM type, with special weight to type III EGMs, which are the distinguishing parameter.

The HT method showed the highest results in terms of a confusion matrix for all EGM types, being considered this way the most reliable method, achieving high performance. On the contrary, the AAT method showed relatively low results in AF type I EGMs, due to the high amplitude variability present in this EGM type of the current database.

#### *4.2. Cycle Length as an Evaluation Parameter*

The main implementation of LAWs detection in real-time mapping devices is related to CL estimation. Consequently, the error in CL estimation was included in the evaluation procedure. As it was shown, mean CL error can be affected by under- or oversensing incidence, and hence, its employment is not recommended. Contrarily, individual CL estimation is an impartial method. Individual CL error was the lowest in the HT method for all EGMs, implying again the high reliability and performance of this technique, while the AAT and ACLI methods have shown poor results, especially for type III EGMs.

#### *4.3. Effect of Signal Length*

Adaptability and proper function regardless of external factors are keys to optimized performance. As signal length varies from study to study, the ability of an algorithm to perform optimally under any length should be taken into consideration for the evaluation of its performance. Additionally, signal length could affect the conclusions of the conducted comparison. For the aforementioned reasons, the individual CL error on each method was assessed over various EGM lengths.

The robustness of the HT algorithm was proved not only by the fact that the error was minimal for any EGM type and length but also because it was not significantly affected by the signal length. This means that this method can perform regardless of the duration of the database used, providing flexibility both to prospective studies and to CA devices. On the other side, the AAT method showed the highest individual error for all EGM types, especially high for EGM type I. Moreover, this algorithm showed unusually high dependence on EGM length for types I and II EGMs, possibly due to high discrepancies between amplitude values in various segments of the same recording.

#### **5. Conclusions**

The HT method is an optimal LAW detector, showing higher performance than the AAT and ACLI techniques in every aspect under investigation. The high accuracy and robustness of the HT method contribute to precise CL length estimation, regardless of the signal length or the EGM fractionation. The high adaptability of this method suggests its implementation on real-time mapping devices aiming to localize CA targets or its adoption by future studies, providing reliable LAW detection results and precise CL estimation.

**Author Contributions:** Conceptualization, D.O., A.V., R.A. and J.J.R.; methodology, D.O., A.V., R.A. and J.J.R.; software, D.O.; validation, D.O., J.M.-A., V.B.-G., R.A. and J.J.R.; resources, J.M.-A., V.B.-G. and J.J.R.; data curation, D.O., J.M.-A. and V.B.-G.; original draft preparation, D.O. and A.V.; review and editing, D.O., A.V., J.M.-A., V.B.-G., R.A. and J.J.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has received financial support from public grants PID2021-123804OB-I00 and PID2021-00X128525-IV0 of the Spanish Government with DOI 10.13039/501100011033 jointly with the European Regional Development Fund (EU), and regional grants SBPLY/17/180501/000411 from Junta de Comunidades de Castilla-La Mancha and AICO/2021/286 from Generalitat Valenciana.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, complying with Law 14/2007, of July 3rd, on Biomedical Research and other Spanish regulations and was approved by the Ethical Review Board of the University Hospital of San Juan (San Juan de Alicante, Alicante, Spain) with protocol code 21/046.

**Informed Consent Statement:** Written informed consent was granted from all the subjects participating in the present research. All acquired data were anonymized before processing.

**Data Availability Statement:** The data supporting reported results and presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors have no association with commercial entities that could be viewed as having an interest in the general area of the submitted manuscript. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


*4.6. Prospects for the Further Development of Ultrasound Methods and New Devices for Milk Analysis*

New promising directions for the development of milk characterization can become microfluidic concepts, phononic crystals concepts, and a combination of acoustic and optical approaches in one device.

Microfluidic technologies may be of interest for separating milk ingredients using acoustophoresis and other methods. Pre-separation of particles based on their differences in size and other physical properties allows for more accurate analysis at a later stage. Up to now, various microfluidic devices have been developed to separate and characterize suspensions and emulsions [112,129]. The acoustofluidic separation can become a powerful instrument for accurate milk analysis.

Improving existing and developing new microfluidic sensor devices involves the use of phononic crystals with liquid-filled cavities [146–148]. Phononic crystals can manipulate and tailor elastic and acoustic wave propagation, can create dynamically resonant metastructures. Phononic crystal designs can be scaled to meet given real object dimensions. Resonance frequency and bandwidth are highly sensitive to geometric and material parameters, specifically speed of sound and sound attenuation. High resolution requires high-Q resonators. Fabrication of phononic crystal and resonant sensor structures with high Q-factor is a challenging but realistic task. Recent studies [146–148] show that the use of phononic crystals can improve the resolution and isolation of responses from liquid resonances. The unique combination of the properties and concepts of phononic crystals, resonant sensors, ultrasonic measurement, and acoustofluidic methods can lead us to a new generation of high-sensitive and high-resolutive liquid sensor devices.

The development of liquid sensor devices that combine the advantages of optical and acoustic methods are carried out both on the basis of traditional methods (like ultrasonic interferometer include optical diffraction [141]) and using promising phononic-photonic crystals approaches [149].