Abstract

Information about potential earthquake sources is a key issue for seismic hazard assessment. This study presents the application of a phenomenological approach based on pattern recognition to determine the possible locations of strong earthquakes in the Bulgarian region. The technique assumes the origin of strong earthquakes in morphostructural nodes formed around the intersections of morphostructural lineaments identified by morphostructural zoning. For the territory of the Bulgaria and neighbouring regions, 178 nodes were defined in this work. The CORA-3 pattern recognition algorithm identified 59 seismogenic nodes analysing a set of geophysical and geological node’s characteristics. The identified seismogenic nodes are capable to generate earthquakes with magnitude equal to or greater than 6 and are located at the boundaries between the largest tectonic domains: Rila, Pirin, and Rhodope orogens; the Serbian-Macedonian massif; and in the Stara Planina belt. The set of characteristic features of seismogenic nodes indicates that the vicinity of potential nodes is characterized by a high contrast of neotectonic movements of the Earth’s crust and the presence of deep heterogeneities in the Earth’s crust. About 40% of the recognized nodes are not associated with any earthquakes, while the rest of the recognized seismogenic nodes are characterized by an area with a radius of 25 km where earthquakes are known to occur. Part of these “non active” seismogenic nodes are close to the historical events with magnitudes higher than 5.5 since the magnitude and location of historical events have large uncertainties. Another part of the seismogenic nodes may slightly change the location due to the uncertainties in morphostructural zonation. Other nodes may indicate unknown historical seismicity or paleoearthquakes. Defined M6+ seismogenic nodes can fill the potential gaps in the recorded seismicity on the territory of Bulgaria, thus to improve the seismic hazard assessment of the studied region.

1. Introduction

Over the centuries, the Bulgarian region has been exposed to a high seismic activity [1, 2]. Precise information about the sources of earthquakes (mainly location and magnitude) significantly improves the results of the seismic hazard assessment, obtained by variety of methods.

Over the past decades, many seismic hazard assessment studies have been carried out for the Bulgarian region [39]. Most of these works implemented the probabilistic seismic assessment approach (PSHA), the reliability of which has come under considerable criticism in recent years [10, 11]. Earlier, the Carpathian-Balkanides mountain belt was studied by this morphostuctural zonation and pattern recognition approaches, and several seismogenic nodes were identified within the Stara Planina region [12, 13]. It involves the hypothesis of the location of strong earthquake epicenters in morphostructural nodes formed at the intersection of morphostructural lineaments limiting morphologically homogeneous blocks of the Earth’s crust [1416]. Specialized morphostructural zonation (MSZ) [15, 17] allows the determination of node position on the Earth’s surface. The employed hypothesis on the triggering of strong earthquakes at the nodes initially was based on observations [14]. There are a number of theoretical works that explain stress concentration around intersections of the faults. Theoretical studies provide the physical explanation of the discussed phenomena. Das and Aki [18] proposed a model for the development of a seismic rupture on a fault, the plane of which is complicated by various kinds of barriers. The barriers are responsible both for stopping the propagation of the rupture along the fault plane and for initiating the seismic rupture, due to the concentration of tectonic stresses in them [19]. Fault zone intersections provide locations for the initiation and healing of ruptures, according to King [20]. Movements along one of the intersecting faults are blocked by another fault; therefore, sufficient stress for the origin of an earthquake can be accumulated at the intersection, as shown by Talwani [21].

The neodeterministic seismic hazard assessment approach (NDSHA) was used in some works for certain regions of Bulgaria, but these studies did not cover the entire territory of Bulgaria and did not use data on potential seismogenic nodes [22]. In particular, the importance of seismogenic nodes in NDSHA applications has been demonstrated in several researches [23, 24]. Long-term studies of seismically active regions of the world have shown a fairly high reliability of the used methodology (morphostructural zonation and pattern recognition). About 86% of postpublication earthquakes occur at seismogenic nodes that were previously recognized as seismically hazardous [25, 26].

Information on possible locations of strong earthquakes obtained using the morphostructural zonation and pattern recognition methodology can be directly incorporated into the NDSHA procedure for seismic hazard assessment, thereby filling possible gaps in known seismicity [27].

In the present study, the potential locations of earthquakes with magnitude equal to or greater than 6 (M6+) are identified in the Bulgarian region, the geographical scheme of which is shown in Figure 1. The magnitude threshold corresponds to the earthquakes reported by the earthquake catalogues in Bulgaria [1, 2] causing damages and casualties.

2. Methodology

The phenomenological approach based on pattern recognition is applied in order to identify potential locations of strong earthquakes. The results of this study provide new information that allow the reduction of the possible uncertainties associated with seismic source identification, which is the first and important element of any seismic hazard assessment, and may improve the performance of future seismic hazard maps for the Bulgarian region.

The employed approach is based on two principal steps. The first one essentially consists of the analyzed object definition—the morphostructural nodes—by morphostructural zonation method (MSZ). The second step uses a pattern recognition algorithm to classify the defined MSZ nodes into two groups: nodes capable to generate earthquakes with magnitude above a defined limit and nodes where seismic events with smaller magnitude may occur only.

MSZ is employed to define a system of areas that are hierarchically ordered. Each area is characterized by homogeneous quaternary topography and uniform tectonic structure. MSZ delineates lineament-and-block structure that includes the following three elements: (1) block areas (blocks) of specific rank (I, II, and III); (2) morphostructural lineaments—block’s boundary zones with different rank; (3) nodes that are intersection of lineaments. MSZ does not take into account information about the seismicity of the region in order not to influence the definition of MSZ objects.

The nodes are formed at sites where block boundaries of different orientation intersect and are characterized by particularly intensive fracturing and contrasting neotectonic movements with a mosaic pattern of structure and topography resulting [17, 28].

The pattern recognition technique assumes that the nodes situated close to the epicenters of one or more strong earthquakes have similar characteristics. This suggestion is used to define seismogenic nodes that have the same features, but any strong seismic event is not identified yet in their vicinity. The pattern recognition separates all nodes into two classes: (1) class D (dangerous nodes)—earthquakes with magnitude may occur in their vicinity; (2) class N (nondangerous nodes)—capable to generate only earthquakes with , where denotes the magnitude threshold of the target earthquakes for a study region. In this work, we fix equal to 6.0 because events of this size regularly occur in the region and expose the real threat for population and economy.

The nodes are characterized by a set of topographical, geological, and geophysical parameters, represented as the scalar components of vectors, which are unique characteristics of each node without any connection with seismicity data in a studied region. The set of all nodes' vectors was used as an input data in pattern recognition algorithm CORA-3 [14, 15].

3. Morphostructural Zoning of the Bulgarian Region

The area of interest is characterized by a complex and heterogeneous structure that has developed during a long tectonic evolution [2931]. Since the Cenozoic, the region has been dominated by an extensional regime [32]. According to Dabovsky et al. [33], the territory of Bulgaria consists of two main tectonic domains: the northern part of the Alpine thrust belt in the Balkans and the Moesian platform adjacent to the north. These major domains are divided into a number of tectonic zones. Three extended tectonic zones are distinguished within the Balkan orogenic system: Balkan, Srednogorie, and Moravian-Rhodope. They are composed of rock complexes of different age and composition [34] and are separated by large regional faults. These main faults are dissected by a dense network of smaller faults, and some of them are classified as active [35, 36]. The thickness of the Bulgarian crust varies in the range of about 30 km near the Black Sea and more than 50 km at the southwest part of the region [37, 38].

The morphostructural map is based on the joint analysis of different maps (topographic, tectonic, and geological), satellite photos and publications on geomorphology, active faults, and geology of the studied region.

MSZ specifies a three-level hierarchy of morphostructures [15, 17]. The lower IIIrd rank unit in the hierarchy is block, an area with a similar value of the quantitative indices of large landforms. Blocks build up the IInd rank units called megablocks. A common trend in variations of elevation and orientation of large landforms characterizes megablocks. Megablocks are united into Ist rank units, mountain countries, characterized by a uniform topography and tectonic structure. Quantitative indices of large landforms such as their height and orientation, as well as the trend of the variations of these indices, are analyzed to delineate blocks and megablocks. A morphostructural lineament is treated as a boundary zone between areal units delineated by MSZ. Morphostructural lineaments are also classified as territorial units; their rank depends on the rank of the unit bounded by the lineament. Taking into account the regional tectonic orientation and the predominant strike of the relief, the boundary zones are defined as longitudinal or as transverse lineaments. The regional strike of the topography and tectonic structure is approximately parallel to the longitudinal lineaments that delimit major landforms and usually follow prominent faults. Transverse lineaments cross major landforms ,and they are discontinuously in the relief, represented by ledges, linear segments of river valleys, and pieces of faults lying in the same strike. Nodes are areas of some extent formed around the intersections of lineaments. Lineaments manifest as zones of a specific width on the Earth’s surface. A real size of the lineament zone can be mapped from fieldwork. According to [17, 39], the width of the Ist rank lineaments, for example, varies from 5 km to 25 km. Additionally, the width of this complex zone may vary along a given lineament. Figure 2 displays the MSZ blocks and lineaments of Ist and IInd ranks along with the name of the units used in this study.

First-rank units and their boundaries (Ist rank lineaments) are identified by MSZ, and usually, the boundaries separate the study region from the neighbouring large tectonic units (Figure 2). In the west, the lineament of the Ist rank 4-34-76-170 separates the study region from the Dinarides and Hellenides. In the south, lineament 4-7-8-13 of the Ist rank separates the Rila-Rhodope crystalline massif from the Aegean Sea basin. In the east, the lineament of the Ist rank 13-98-198 separates the continental morphostructures from the deep Black Sea basin. This boundary is traced along the continental slope of the Black Sea. To the north, the MSZ map is bounded by lineament 175-191-198 of the IInd rank, passing along the Danube River. According to the MSZ rules [15, 16, 40], the Bulgarian region was divided into six units of the Ist rank (so-called mountain countries): southern part of the Moesian platform, Stara Planina, Srednogorie, Serbia-Macedonia massif, Rila-Rhodope massif, and Thrace basin (shown in Figure 2). A detailed description of several mountain countries is presented in the following text.

The Moesian platform is composed of relatively undeformed Mesozoic rocks. It is divided into two megablocks by a transverse lineament of the IInd rank 146-191. The megablocks differ in relief and drainage patterns. The topographic elevation is higher to the east of the lineament, and the predominant river orientation becomes NW-SE, while the western part of the Moesian platform is dominated by a NE-SW river orientation.

Stara Planina, represented by an extended linear orogen, is part of the Carpathian-Balkan arc. To the north, it borders the Moesian platform. On the map of the MSZ, this boundary is represented by a lineament of the Ist rank corresponding to the Fore-Balkan fault 99-146-161-175. To the south, a lineament of the Ist rank separates the relatively high Balkanides from the Srednogorie, which is characterized by lower elevations. The zone of this lineament corresponds to the sub-Balkan fault 98-134-172. The eastern boundary of the Stara Planina is a transverse lineament of the Ist rank, traced in accordance with the cryptostructure of the Black Sea 98-99. Transverse lineaments of the IInd rank have both NW-SE and NE-SW orientations. Lineaments of the NE-SW orientation are consistent with transcurrent faults identified by geological and geophysical methods [30, 41]. Based on the analysis of satellite data, a system of NW-SE lineaments of the Balkan region, including the Balkanides, is identified in [29].

Srednogorie is located between the Stara Planina mountain region and the Rila-Rhodopean massif. To the north, Ist rank lineament 98-134-172 separates Srednogorie from Stara Planina; the lineament corresponds to the sub-Balkan fault [42]. To the south, Ist rank lineaments 80-88-20 and 45-14 are partly consistent with the Maritsa strike-slip fault zone [43, 44]. In the relief of the Srednogorie, an elevated western and central part and a flat eastern part stand out. In accordance with the change in the height of the relief, three megablocks were identified which are divided by transverse lineaments 86-136, 55-107, and 52-97 of the IInd rank into three megablocks: western, central, and eastern. This division is consistent with the major structural units identified in the Srednogorie by Dabovsky et al. [33].

The rest of the mountain countries have much more complex structure.

Lineaments of the IIIrd rank divide megablocks into blocks, territorial units of the IIIrd rank. The blocks are characterized by a similar height and orientation within the same landform. Lineaments of the IIIrd rank control local changes in the height and orientation of large landforms along their strike.

The compiled MSZ map (Figure 3) presents a network of lineaments and a hierarchical system of blocks within the Bulgarian region. Relative displacements and interaction of these blocks may cause shallow seismicity in the region.

4. Recognition of Seismogenic Nodes M6+ in the Bulgarian Region

In total, 178 intersections of lineaments are identified in the studied region. Since the morphostructural map is compiled on the basis of cartographic sources without field research, natural boundaries of the nodes had not been delimited. Previously, in the Caucasus region, the nodes were mapped during a fieldwork, and the node’s dimensions were evaluated from 20 to 60 km in radius [28, 39]. Thus, the areas within the radius of 20-30 km around the points of lineament intersections can be regarded as nodes. Within the present study, nodes are assumed to be circles of a 25 km radius around the points of lineament intersections. Such valuation proves to be in agreement with the size of earthquake source for M6+ [45, 46].

4.1. Learning Sets Selection

Seismicity data used to assemble the training set is presented in Table 1. The main events with M6+ have been selected from the Bulgarian earthquake catalogues [1, 2]. Application of the CORA-3 algorithm requires a training set of objects (vectors), for which we assume a priori the class (D or N) they belong to. This training set is formed by three subsets: , , and . In the case of the Bulgarian region, the training objects of subset are all nodes hosting earthquakes with that occurred after 1900 (see Table 1). The training objects of subset include the nodes that are not related to any events from Table 1. The subset includes nodes marked by earthquakes M6+ that happened before 1900. CORA-3 uses only nodes from the subsets and for the selection of the characteristic traits. Subset is processed at the recognition stage, when each node is assigned either to D or to N class.

As a result, the subset consists of 16 objects, the subset includes 142 objects, and the subset is formed by 20 objects (Figure 4). It is important to emphasize that is not a “pure” training set because some of its members may belong to class D, i.e., earthquakes with may occur near some of these objects but are not known at present, since the period of observations is short.

4.2. Parametrization of Recognition Objects

Objects of recognition are presented by the morphostructural nodes. MSZ delineated 178 nodes. Each object is described by a vector, the components of which are values of parameters measured uniformly for all 178 nodes. The goal of the pattern recognition is to divide the vectors (and their corresponding nodes) into classes D and N.

The parameters of the nodes, used in this work, were estimated in a circular area with a radius of 25 km centred at each node. In total, 18 parameters have been considered in this work, and they are presented in Table 2. Morphometric parameters have been measured from topographic maps at the scale of 1 : 500 000. Gravity and magnetic characteristics were defined using the GeoMapApp database (http://www.geomapapp.org).

After measuring the values of these 18 parameters for all nodes, the corresponding vectors with 18 components are considered as recognition objects or recognition patterns.

These parameters are based on the geomorphic and morphometric information, the characteristics of the block-and-lineament geometry, and the gravity and magnetic data. Morphometric parameters characterize the contrast and intensity of tectonic deformations during the period of the relief building. Parameters of the block-and-lineament geometry reflect the complexity and fragmentation of the media in the node vicinities. Gravity and magnetic data characterize the deep heterogeneity of the crust. In principle, other relevant information characterizing the specific features of seismic areas may be considered for recognition. The only necessary precondition is that the value of each parameter should be defined with the same accuracy for each recognition pattern (i.e., the node) within the study region. Most of the parameters used to recognize D nodes in the Bulgaria region were used in other regions previously studied with this approach [1217, 24, 25, 28, 39, 47]. A special feature of this work is the use of the most complete set of parameters to date as compared to the majority of the previous studies. This applies primarily to parameters based on gravimetric and magnetic data.

The CORA-3 algorithm works in a binary vector space. Therefore, prior to running the algorithm, the objects of recognition that are vectors with real components (actual parameter’s value determined from maps or digital databases) should be converted in vectors with binary components. For this purpose, the ranges of the parameter’s value are divided into by one or two thresholds that defined two or three disjoint ranges, correspondingly. Thresholds are specified by the discretization procedure described in details in several papers [14, 15, 28]. If a range of a parameter is divided into two parts by one threshold, then only two gradations (“small” and “large”) are used after the discretization instead of the actual parameter’s value. In this case, “small” values are coded as 1, and “large” ones are coded as 0. If the range is divided into three parts by two thresholds, then the number of intervals is three: “small” (11), “medium” (01), and “large” (00).

The discretization causes loss of information but makes the results of recognition more stable to the changes and variations of the parameters data. After the discretization, the parameter values are converted into components of binary vectors, where each vector has 18 binary components corresponding to each object of recognition.

4.3. Recognition of the Seismogenic Nodes

There are 39 known main earthquakes M6+ (see Table 1) for the period from the Ist century to the end of the XXth century on the studied territory [1, 2]. Objects of recognition are associated with the morphostructural nodes. The total number of nodes is 199, but during processing, a few of them were omitted, resulting in a final utilization of 178 nodes (Figure 3). Each of the 178 nodes is described by a parameter vector. The problem to be solved is to separate the vectors (and the respective nodes) into classes D and N. The CORA-3 algorithm is applied to 178 vectors (as objects of recognition) to classify them into classes D and N. The decision rule for the classification is composed by 12 characteristic traits of class D (D-traits) and 11 characteristic traits of class N (N-traits) that were selected by CORA-3. They are given in Table 3 in the form of inequalities on the values of the parameters in accordance to the discretization thresholds (Table 2). The classification of all nodes into D and N is made as follows. For each node, the algorithm calculates the number of the characteristic traits for class D, the number of those for class N, and the difference . Class D includes the nodes for which , while class N consists of the nodes for which . The value of is specified as 1.

Twelve D-traits and eleven N-traits have been selected by the algorithm. In total, CORA-3 algorithm assigned 59 nodes out of 178 delineated nodes to seismogenic class D. These nodes are shown by gray circles in Figure 5. The list of the recognized D nodes is presented in Table 4.

4.4. Control Experiments

The results obtained by pattern recognition for the territory of the Bulgarian region (presented in Figure 5) are called “main results” in the following control experiments. A set of control tests is executed to evaluate indirectly the reliability of the classification obtained by the recognition algorithm [14, 15, 28].

In the test called “Seismic Future,” the subsets and include all objects recognized in the main variant as D and N, respectively. Classification of 10 objects (about 6% of the total number of objects) is changed in this test compared with the main variant.

The test “Objects exclusion” allows to estimate contributions of the individual parameters listed in Table 2 to the main result. For this purpose, each of the eighteen parameters governing the main variant was eliminated from consideration one at a time, and classification was repeated by the CORA-3 algorithm. In this test, 12 objects (7% of the total number of objects) changed their classification.

In the test “Function exclusion,” the parameters composing the decision rule (Table 3) were also excluded in turn from the recognition procedure. About 15 objects (around 8%) changed their classification with respect to the main result.

According to Gvishiani et al. [28], such insignificant changes in the main result indicate its sufficient stability.

5. Discussion

Most of the identified earthquake-prone nodes are located on the lineaments of the highest (Ist and IInd) ranks, which correspond to the prominent faults and divide the largest blocks of the Earth’s crust in the Bulgarian region. The majority of D nodes are located in the mountainous areas, i.e., in Rila, Pirin, and Rhodope massifs, as well as in the Stara Planina Mountain range (see Figure 5). In general, the spatial distribution of seismogenic nodes corresponds to the dispersion of the observed seismicity in the Bulgarian region. Most of the D nodes are identified in the western part of Bulgaria, where the observed seismicity is concentrated. The obtained results also confirm the high seismic potential of the Shabla region (northeast Bulgaria). Nodes 101 and 100 are also identified as seismogenic in the Shabla area, including nodes 99 and 149 located in the Black Sea (see Figure 5). Potential nodes (46 and 15) were identified on the boundary between the Strandzha area and the Thrace basin where M6+ earthquakes were not observed up to now; and in the node’s vicinity, there are no active faults according the European Database of Seismogenic Faults (EDSF) [48]. The possible reason could be an inaccuracy in the MSZ or unknown paleoearthquakes or unknown faults.

In the most northwestern part of the region, node 170, in which two historical events of 1739 and 1895 took place, is not recognized as D type. This node is located in the Serbo-Macedonian massif, which differs in relief and structure from those in the Bulgarian region. This node is recognized as capable of generating earthquakes of M6+ in the similar study for the Dinarides [47]. Figure 6 displays the recognized seismogenic nodes including node 170 as seismogenic for M6+.

In agreement with characteristic features defined by CORA-3, D nodes are characterized by “large” values of the relief energy and relief gradient in combination with “large” values of gradients of gravitational and magnetic anomalies (Table 3). The important role of gravity and magnetic anomalies in the spatial distribution of seismicity in Bulgaria was shown by Trifonova et al. [49]. In general, the set of characteristic features of seismogenic nodes established by recognition indicates the high contrast of neotectonic movements and the presence of deep inhomogeneities in the Earth’s crust [31] in the vicinity of D nodes.

The statistics of the obtained results shows that 33% of the nodes on the Bulgarian territory are seismogenic for events M6+ (59 from totally 178 nodes). In the 25 km radius of 14 out of 59 nodes, there is at least one documented event after 1900. Close to 20 nodes out of 59, there is at least one known event before 1900. Earthquakes of magnitude 5.5 to 5.9 before 1900 were observed in a 25 km vicinity of 9 nodes out of 59 nodes since the location and magnitude of historical events have large uncertainties (some events from [1] are revised in Glavcheva and Matova [50]). All 59 nodes of class D are in areas with geological faults according to a geological map of Bulgaria with a scale of 1 : 100 000 [51]. Fifty-six of the seismogenic nodes are located very close to active faults from the EDSF database [48] (Figure 6).

The performed recognition technique identifies a number of D nodes where events M6+ have not been documented up to present. This is a significant amount of additional information in order to improve seismogenic source models and contribute to long-term seismic hazard assessment in the Bulgarian region. The information on the possible locations of strong earthquakes provided by this work can be directly incorporated in the neodeterministic procedure for seismic hazard assessment, thus, filling in possible gaps in known seismicity in the Bulgarian territory.

There are several events with M7+ in the Bulgarian region, but the number of nodes to which they are associated proves to be insufficient for a proper training of the algorithm. It is possible to solve the recognition problem for M7+ events studying larger territory (for example, Bulgaria studied in the paper of Gorshkov et al. [52] and Greece studied by Novikova and Gorshkov [53]), in similarity to the work of Novikova and Gorshkov [54], where three neighbouring regions (Caucasus, Kopet Dag, and Alborz) were jointly considered to recognize dangerous nodes for M7+.

6. Conclusions

Using the morphostructural zonation technique, we define 178 nodes for the territory of Bulgaria and neighbouring regions. A set of geophysical and geological data are assigned to each node to characterize its properties. The CORA-3 pattern recognition algorithm uses a subset of nodes to “train” identification of seismogenic nodes capable to generate earthquakes with magnitudes equal to or greater than 6. Most of the recognized seismogenic nodes (59 nodes) are located at the boundaries between the largest tectonic domains: the border between Rila, Pirin, and Rhodope orogens; the Serbo-Macedonian massif; and in the Stara Planina belt. The seismogenic nodes in the region are characterized by the high contrast of neotectonic movements of the Earth’s crust and the presence of deep heterogeneities in the Earth’s crust. The area of about 40% of the recognized nodes has not yet observed any earthquake M6+. Several conclusions on these still “silent” nodes can be made: they are close to the historical events with magnitudes between 5.5 and 6.0 since the magnitude values have a large uncertainty; they may slightly change the location due to the uncertainties in morphostructural zonation; they may be a sign of past earthquake activity that has not been documented because the area is not populated.

The upper threshold of possible earthquakes in the region cannot be resolved by the applied methods. This problem can be solved, for example, on the basis of seismicity modeling based on a block-and-fault model, as was done for other regions of the world [55, 56]. The preliminary results of block-and-fault seismicity modeling for the territory of Bulgaria are shown in Dimitrova et al. [57], but the work on improving the results is in progress.

The information on the possible locations of strong earthquakes obtained within the frame of this work can be directly incorporated in the neodeterministic seismic hazard assessment, complementing the documented seismicity on the Bulgarian territory with a number of probable strong earthquake locations (past unknown and future events).

Data Availability

The data used in this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Authors’ Contributions

Alexander Gorshkov was responsible for the conceptualization, methodology, software, formal analysis, and writing of the original draft preparation. Olga Novikova was responsible for the software, formal analysis, writing of the original draft preparation, and project administration. Sonya Dimitrova was responsible for the investigation, resources, data curation, and writing of the original draft preparation. Lyuba Dimova was responsible for the resources, data curation, writing of the review, and editing. Reneta Raykova was responsible for the conceptualization, formal analysis, investigation, writing of the review, editing, and funding acquisition. All authors read and approved the final manuscript.

Acknowledgments

This study is funded by the Russian Foundation of Basic Research (RFBR) (research project no. 20-55-18008) (partly) and by the Bulgarian National Science Fund (research project no. KP-06-Russia-29/16.12.2020).