*Article* **Bioinspired Electropun Fibrous Materials Based on Poly-3-Hydroxybutyrate and Hemin: Preparation, Physicochemical Properties, and Weathering**

**Polina M. Tyubaeva 1 , Ivetta A. Varyan 1,2, Anna K. Zykova 1 , Alena Yu. Yarysheva 3 , Pavel V. Ivchenko 3,\* , Anatoly A. Olkhov 1,2 and Olga V. Arzhakova 3**

**Abstract:** The development of innovative fibrous materials with valuable multifunctional properties based on biodegradable polymers and modifying additives presents a challenging direction for modern materials science and environmental safety. In this work, high-performance composite fibrous materials based on semicrystalline biodegradable poly-3-hydroxybutyrate (PHB) and natural ironcontaining porphyrin, hemin (*Hmi*) were prepared by electrospinning. The addition of *Hmi* to the feed PHB mixture (at concentrations above 3 wt.%) is shown to facilitate the electrospinning process and improve the quality of the electrospun PHB/*Hmi* materials: the fibers become uniform, their average diameter decreases down to 1.77 µm, and porosity increases to 94%. Structural morphology, phase composition, and physicochemical properties of the *Hmi*/PHB fibrous materials were studied by diverse physicochemical methods, including electronic paramagnetic resonance, optical microscopy, scanning electron microscopy, and energy-dispersive X-ray spectroscopy, elemental analysis, differential scanning calorimetry, Fourier-transformed infrared spectroscopy, mechanical analysis, etc. The proposed nonwoven *Hmi*/PHB composites with high porosity, good mechanical properties, and retarded biodegradation due to high antibacterial potential can be used as high-performance and robust materials for biomedical applications, including breathable materials for wound disinfection and accelerated healing, scaffolds for regenerative medicine and tissue engineering.

**Keywords:** biodegradable polymers; poly(3-hydroxybutyrate); hemin; electrospinning; outdoor weathering; new technologies

## **1. Introduction**

At the present time, the mainstream direction of modern materials science and technology is oriented toward the development of innovative and ecologically friendly materials based on natural components and biodegradable polymers with task-oriented multifunctional properties presents. In this connection, the family of biodegradable polyhydroxyalkanoates, including poly(hydroxybutyrate) (PHB) is of special interest due to their potential as biomedical materials for therapeutic applications. PHB seems to be a promising candidate for the development of composite materials for biomedical applications as this biopolymer is characterized by high melting temperature, high degree of crystallinity, and low permeability to oxygen, water, and carbon dioxide [1]. This biodegradable polymer can be synthesized from renewable natural sources and is known to be biocompatible with the human body [2].

The introduction of various additives to PHB allows the preparation of the composite PHB-based materials with desired functional properties (nanoparticles [3], carbon

**Citation:** Tyubaeva, P.M.; Varyan, I.A.; Zykova, A.K.; Yarysheva, A.Y.; Ivchenko, P.V.; Olkhov, A.A.; Arzhakova, O.V. Bioinspired Electropun Fibrous Materials Based on Poly-3-Hydroxybutyrate and Hemin: Preparation, Physicochemical Properties, and Weathering. *Polymers* **2022**, *14*, 4878. https://doi.org/ 10.3390/polym14224878

Academic Editor: Piotr Bulak

Received: 23 September 2022 Accepted: 9 November 2022 Published: 12 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**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/).

nanotubes [4], catalysts, and enzymes [5]). Nowadays, special interest is focused on the molecular complexes of porphyrins and related materials for their applied use as superconductors [6–8], in analytical chemistry [9,10], photonics [11,12], gene therapy [13,14], biomedicine [15–17]. Progress in the synthesis and chemistry of tetrapyrroles has triggered the ever-growing attention to tetrapyrrole complexes as an effective modifying additive for polymeric materials, including numerous analogs of natural systems such as porphyrins, texaphyrins, chlorins, corroles, and others, applied for biomedical purposes [18–22].

Hemin is the iron-containing porphyrin, which is endogenously produced in the human body, for example, during the turnover of old red blood cells. The unique properties and benefits of hemin have been widely discussed in the literature [23–32]. Evidently, hemin is biocompatible with living organisms [23] and can be used for the stabilization of protein molecules on polymer carriers, and hemin-containing materials can be successfully used for diverse biomedical applications [24], including, for example, the development of safe containers for drug delivery systems [25]. Moreover, peptides can be attached to the surface of hemin-containing materials [26,27]. Hemin shows a well-pronounced catalytic activity [28] and can be effectively used in anticoagulants [29]. Hemin is characterized by high thermal stability, thus lifting the challenging problems related to the processing of polymers and the preparation of sustainable polymer nanocomposite materials [30]. In addition, hemin also shows high antimicrobial activity against Gram-positive *S. aureus* [31], which is a significant advantage in the creation of biomedical materials. Finally, the benefits of hemin are concerned with its natural origin biocompatibility with living organisms, and simple and facile procedure of synthesis [32]. Hence, hemin was selected as an efficient additive for the modification of the PHB materials and preparation of new *Hmi*/PHB materials with task-oriented functional properties.

Electrospinning (ES) is credited as a reliable and efficient method for the preparation of high-performance nonwoven materials based on a wide range of polymers [33]. This method provides the formation of an interconnected network of ultrathin fibers by drawing a jet of the polymer solution under the action of physical forces. Electrospinning offers powerful tools for the preparation of diverse composite materials when a functional additive is dissolved in the polymer feed solution [34–37]: poly(3-hydroxybutyrate-co-3-hydroxyvalerate)/polylactide [38], poly(butylene adipate succinate)/polylactide [39], poly(butylene adipate terephthalate)/polylactide [40], poly(ε-caprolactone)/poly (3-hydroxybutyrate) [41]. This ES approach allows one to control the content of additives in the composite, phase composition, and distribution of the additive within the fibrous materials.

This work addresses the development of a facile and sustainable method for the preparation of high-performance innovative biocompatible *Hmi*/PHB materials by electrospinning, characterization of their morphology, phase composition, and mechanical properties by diverse physicochemical methods (electronic paramagnetic resonance (EPR), optical microscopy, scanning electron microscopy (SEM) and energy-dispersive X-ray spectroscopy (EDX) elemental analysis, differential scanning calorimetry (DSC), Fouriertransformed infrared (FTIR) spectroscopy, mechanical analysis, etc.), investigation of their stability upon outdoor weathering tests, and description of the potential applied benefits of the resultant *Hmi*/PHB materials for diverse applications.

#### **2. Materials and Methods**

#### *2.1. Materials*

In this work, semicrystalline biodegradable polymer, poly(3-hydroxybutyrate) (PHB) (trade mark 16F series, BIOMER, Frankfurt, Germany) was used; molecular weight— 206 kDa, density—1.5 g/cm3 and degree of crystallinity—59% (Figure 1A). Hemin (a tetrapyrrole complex from the class of natural porphyrins) was used as a modifying additive (Figure 1B). Hemin with the coordination complex of iron (oxidation state III) was obtained from the bovine (Aldrich Sigma, Saint Louis, MO, USA).

obtained from the bovine (Aldrich Sigma, Saint Louis, USA).

**Figure 1.** Structural formulae of PHB (**A**) and hemin (**B**). **Figure 1.** Structural formulae of PHB (**A**) and hemin (**B**).

#### *2.2. Preparation of Electrospun Materials*

*2.2. Preparation of Electrospun Materials* In this work, polymer nanofibrous materials were prepared by electrospinning (ES) technique using a single-capillary laboratory unit with a capillary diameter of 0.1 mm (EFV-1, Moscow, Russia) (Supplementary files, Figure S1). As a result of electrospinning, polymer fibers are cured due to the complete evaporation of the solvent. For the preparation of feed solutions, finely dispersed PHB powder was dissolved in chloroform at a temperature of 60 °C. Hemin was dissolved in N,N-dimethylformamide at 25 °C and homogenized with the PHB solution. The content of PHB in the solution was 7 wt.%, and the content of hemin was 1, 3, and 5 wt.% of the PHB. The conditions of the ES process were the following: voltage 17–20 kV, the distance between the electrodes—190–200 mm, the gas pressure on the solution—10–15 kg(f) cm−2 . The electrical conductivity of the feed solution was 10–14 μS cm−1 and the viscosity of the hemin/PHB solution was 1.4–1.9 Pa s (viscosity of the 7 wt.% PHB solution in chloroform was 1 Pa s) and the flow rate of the solution was 0.15–0.20 cm3 min−1 (flow rate of the 7 wt.% PHB solution in chloroform was In this work, polymer nanofibrous materials were prepared by electrospinning (ES) technique using a single-capillary laboratory unit with a capillary diameter of 0.1 mm (EFV-1, Moscow, Russia) (Supplementary Files, Figure S1). As a result of electrospinning, polymer fibers are cured due to the complete evaporation of the solvent. For the preparation of feed solutions, finely dispersed PHB powder was dissolved in chloroform at a temperature of 60 ◦C. Hemin was dissolved in N,N-dimethylformamide at 25 ◦C and homogenized with the PHB solution. The content of PHB in the solution was 7 wt.%, and the content of hemin was 1, 3, and 5 wt.% of the PHB. The conditions of the ES process were the following: voltage 17–20 kV, the distance between the electrodes—190–200 mm, the gas pressure on the solution—10–15 kg(f) cm−2 . The electrical conductivity of the feed solution was 10–14 µS cm−1 and the viscosity of the hemin/PHB solution was 1.4–1.9 Pa s (viscosity of the 7 wt.% PHB solution in chloroform was 1 Pa s) and the flow rate of the solution was 0.15–0.20 cm3 min−1 (flow rate of the 7 wt.% PHB solution in chloroform was 0.45 cm3 min−1 ).

kDa, density—1.5 g/cm3 and degree of crystallinity—59% (Figure 1A). Hemin (a tetrapyrrole complex from the class of natural porphyrins) was used as a modifying additive (Figure 1B). Hemin with the coordination complex of iron (oxidation state III) was

#### 0.45 cm3 min−1). *2.3. Methods*

#### *2.3. Methods* 2.3.1. Optical Microscopy

2.3.1. Optical Microscopy Primary data on topography and the general appearance of the samples were col-Primary data on topography and the general appearance of the samples were collected using an Olympus BX43 optical microscope (Tokyo, Japan) in transmission mode.

#### lected using an Olympus BX43 optical microscope (Tokyo, Japan) in transmission mode. 2.3.2. Scanning Electron Microscopy (SEM)

2.3.2. Scanning Electron Microscopy (SEM) The morphology of the samples was examined by SEM observations using a Tescan The morphology of the samples was examined by SEM observations using a Tescan VEGA3 microscope (Wurttemberg, Czech Republic). Prior to SEM observations, the surface of the test samples was decorated with a thin platinum layer.

VEGA3 microscope (Wurttemberg, Czech Republic). Prior to SEM observations, the sur-

#### face of the test samples was decorated with a thin platinum layer. 2.3.3. Estimation of Structural Parameters

2.3.3. Estimation of Structural Parameters The morphology of the fibrous materials was studied using the Olympus Stream The morphology of the fibrous materials was studied using the Olympus Stream Basic software (Tokyo, Japan). The average diameter of individual fibers was estimated for each fiber on the z-stack at five different points. The density δ of the samples reads as:

$$\delta = \frac{m}{l \times B \times b} \tag{1}$$

(1)

*m* = where *m* is the weight of the sample; *l* is the length; *B* is the width; *b* is the thickness.

*l B b* where *m* is the weight of the sample; *l* is the length; *B* is the width; *b* is the thickness. The porosity of the fibrous electrospun materials was estimated as the ratio between the volume of pores and the overall volume of the sample:

$$\mathcal{W} = (1 - \frac{V\_f}{V\_s}) \times 100,\,\% \tag{2}$$

where *W* is the porosity of the nonwoven sample, *Vf* and *Vs* stand for the volume of polymer as fibers in the nonwoven sample and the nonwoven sample, respectively.

#### 2.3.4. Energy-Dispersive X-ray Spectroscopy (EDX)

The elemental composition of the samples (oxygen, chlorine, and iron) was studied by the method of energy-dispersive X-ray spectroscopy (EDX) using a Tescan VEGA3 instrument (Wurttemberg, Czech Republic).

#### 2.3.5. FTIR Spectroscopy

The chemical composition of initial, modified, and weathered materials were studied using an FTIR Lumos BRUKER (Berlin, Germany) spectrometer at a temperature of (22 ± 2) ◦C in the range of wavenumbers of 4000 ≤ ν ≤ 600 cm−6 in the mode of frustrated total internal reflection (FTIR) on the diamond crystal [42]. The intensity of several absorption peaks related to polymer degradation (functional groups) was estimated [43,44].

#### 2.3.6. Differential Scanning Calorimetry (DSC)

The thermal properties of the test samples were studied using a Netzsch 214 Polyma thermal analyzer (Selb, Germany); the heating rate was 10 ◦C/min. The samples were heated from 20 ◦C to 220 ◦C and then cooled down to 20 ◦C. The heat of fusion was calculated using the NETZSCH Proteus software (Selb, Germany). The weight of the samples was ~7 mg. The data were analyzed according to the standard procedure [45]. The degree of crystallinity of the samples *χ* reads as:

$$
\chi = \Delta H\_m / \Delta H\_{m,100\%} \times 100\% \tag{3}
$$

where ∆*Hm* is the heat of fusion; ∆*Hm,*100% is the heat of fusion of the ideal crystal of PHB (146 J/g [46]).

#### 2.3.7. Electronic Paramagnetic Resonance (EPR)

The EPR spectra were collected using an EPR-V automated spectrometer (Moscow, Russia). Microwave power in the resonator cavity was below 7 mW. Stable nitroxyl radical (2,2,6,6-tetramethyl piperidine-1-oxyl (TEMPO) radical moiety) was used as a probe. The TEMPO radical was introduced into the polymer samples from the gas phase at 60 ◦C. The content of the TEMPO radicals in the polymer samples was estimated using Bruker WinEPR and SimFonia software (Berlin, Germany) (CCl4 with the radical concentration below 10−3 mol L−1 as a reference). Correlation times of the probe rotation were estimated from the corresponding EPR spectra.

The EPR spectra of the TEMPO spin probe in the slow-motion region (τ > 10−10 s) were analyzed in the terms of the model of isotropic Brownian rotation using the NLSL software [47]. Correlation times of the probe rotation τ in the fast-motion region (5 × 10−11 < τ < 10−9 s) were estimated from the corresponding EPR spectra as:

$$
\pi = 6.65 \times 10^{-10} (\sqrt{I\_+/I\_-} - 1) \times \Delta H\_+ \tag{4}
$$

where ∆*H+* is the half-width of the spectrum component in the weak field, and *I+/I*− is the ratio of the intensities of the components in weak and strong fields, respectively [48].

#### 2.3.8. Mechanical Tests

Mechanical properties (tensile strength, elongation at break) of the test samples were estimated using a Devotrans DVT GP UG 5 tensile compression testing machine (Istanbul, Turkey) according to ASTM D5035-11.

## 2.3.9. Weathering Experiments upon Outdoor Exposure

The weathering experiments were performed according to the ASTM G7/G7M-21 ("Standard Practice for Natural Weathering of Materials"). The samples were subjected to

outdoor exposure according to standard protocol D1435-13 ("Standard Practice for Outdoor Weathering of Plastics"). According to the protocol, the wooden weathering racks were tilted at an angle of 90◦ with respect to the horizontal line. door Weathering of Plastics"). According to the protocol, the wooden weathering racks were tilted at an angle of 90° with respect to the horizontal line. The dimensions of the samples were 40 mm × 40 mm × 0.3 mm (length, width, and

Mechanical properties (tensile strength, elongation at break) of the test samples were estimated using a Devotrans DVT GP UG 5 tensile compression testing machine (Istanbul,

The weathering experiments were performed according to the ASTM G7/G7M-21 ("Standard Practice for Natural Weathering of Materials"). The samples were subjected to outdoor exposure according to standard protocol D1435-13 ("Standard Practice for Out-

*Polymers* **2022**, *14*, x FOR PEER REVIEW 5 of 19

2.3.9. Weathering Experiments upon Outdoor Exposure

Turkey) according to ASTM D5035-11.

The dimensions of the samples were 40 mm × 40 mm × 0.3 mm (length, width, and thickness, respectively). Similar samples were placed as blank samples in a dark and dry place. The specimens were fastened with white cotton threads; one side of the specimens was fastened with adhesive tape. The duration of each weathering experiment was 3 months (June–August; temperature average 21 ◦C; humidity average 75.7%; geographical location: 55◦ 450 north latitude; 37◦ 370 east longitude). thickness, respectively). Similar samples were placed as blank samples in a dark and dry place. The specimens were fastened with white cotton threads; one side of the specimens was fastened with adhesive tape. The duration of each weathering experiment was 3 months (June–August; temperature average 21 °C; humidity average 75.7%; geographical location: 55° 45' north latitude; 37° 37' east longitude).

#### 2.3.10. Statistical Analysis 2.3.10. Statistical Analysis

All measurements were performed, at least, three times for each sample. Mechanical tests, weathering experiments, and measurements of structural characteristics were performed for 8–10 samples. A mean value and standard deviation were estimated using Data Analysis by the Origin Pro program (OriginLab Corporation, Northampton, MA, USA). All measurements were performed, at least, three times for each sample. Mechanical tests, weathering experiments, and measurements of structural characteristics were performed for 8–10 samples. A mean value and standard deviation were estimated using Data Analysis by the Origin Pro program (OriginLab Corporation, United States).

#### **3. Results 3. Results**

Electrospinning is a well-known and reliable procedure for the development of a fibrous structure [49,50]. This method makes it possible to control the supramolecular structure and morphology of individual fibers as well as the morphology of the system as a whole [51]. The introduction of modifying additives even at low concentrations (below 5 wt.%) provides a marked effect on the formation of the structure of the resultant materials at all structural levels. Electrospinning is a well-known and reliable procedure for the development of a fibrous structure [49,50]. This method makes it possible to control the supramolecular structure and morphology of individual fibers as well as the morphology of the system as a whole [51]. The introduction of modifying additives even at low concentrations (below 5 wt.%) provides a marked effect on the formation of the structure of the resultant materials at all structural levels.

#### *3.1. Morphology Structure and Properties of Hemin/PHB Fibrous Materials 3.1. Morphology Structure and Properties of Hemin/PHB Fibrous Materials*

Nonwoven materials prepared by electrospinning are characterized by well-developed morphology which is described by several key parameters, including density, the average diameter of individual fibers, porosity, and pore dimensions. Figure 2 shows the optical micrographs illustrating the visual appearance of the *Hmi*/PHB samples with different content of *Hmi*. Nonwoven materials prepared by electrospinning are characterized by well-developed morphology which is described by several key parameters, including density, the average diameter of individual fibers, porosity, and pore dimensions. Figure 2 shows the optical micrographs illustrating the visual appearance of the *Hmi*/PHB samples with different content of *Hmi*.

**Figure 2.** The micrographs of *Hmi*/PHB fibrous nanocomposites with different content of hemin: 1 wt.% (**A**), 3 wt.% (**B**), and 5 wt.% (**C**). **Figure 2.** The micrographs of *Hmi*/PHB fibrous nanocomposites with different content of hemin: 1 wt.% (**A**), 3 wt.% (**B**), and 5 wt.% (**C**).

As follows from Figure 2, the *Hmi*/PHB fibrous samples contain big-sized dark inclusions which are accommodated on the surface of individual fibers of the nonwoven sample (Figure 2A, B). The dimensions of the inclusions are equal to 4–32 µm and 0.7–17 µm As follows from Figure 2, the *Hmi*/PHB fibrous samples contain big-sized dark inclusions which are accommodated on the surface of individual fibers of the nonwoven sample (Figure 2A,B). The dimensions of the inclusions are equal to 4–32 µm and 0.7–17 µm for *Hmi*/PHB containing 1 and 3 wt.% of *Hmi*, respectively. Noteworthy is that the *Hmi*/PHB samples containing 5 wt.% *Hmi* are seen to be uniform and free of big-sized inclusions (Figure 2C). This fact indirectly suggests that electrospinning from the feed solutions containing 3–5 wt.% of *Hmi* is enhanced and provides the preparation of more uniform nonwoven *Hmi*/PHB materials. Analysis of the micrographs in Figure 2 allows estimation of the average diameter of individual fibers of the *Hmi*/PHB nonwoven materials. With increasing the content of *Hmi* to 5 wt.%, the average diameter of the fibers decreases from 3.5 µm for the neat PHB down to 1.77 µm for the *Hmi*/PHB containing 5 wt.% of *Hmi*. Noteworthy is that, as the content of *Hmi* in the *Hmi*/PHB composites increases, porosity of nonwoven materials increases from 80% (for neat PHB) to 94%. Hence, this evidence allows one to conclude that electrospinning of the feed solutions with the maximum content of *Hmi* leads to the preparation of nonwoven *Hmi*/PHB materials with uniform structure and high porosity. Structural characteristics of the nonwoven *Hmi*/PHB materials are listed in Table 1. The histogram illustrating the size distribution of PHB and *Hmi*/PHB nonwoven composites is given in the Supplementary Files, Figure S2.

**Table 1.** Structural characteristics of PHB and *Hmi*/PHB nonwoven composites.

Hence, the presence of 5 wt.% *Hmi* in the *Hmi*/PHB feed mixture is associated with a lower density of the feed solution and facilitates the ES process due to better-organized movement through the jet. As a result, pore size decreases, and porosity increases. This is due to the fact that the diameter of the PHB fibers is higher than that of the Hemin-modified fibers. Thicker PHB fibers are also characterized by the presence of large defects and their diameter exceeds the average diameter of the fibers by 3–10 times. As a result, upon ES, thicker fibers are located at a higher distance from each other, forming large pores, and this effect is the consequence of low electrical conductivity of the spinning solution. However, as Hemin is added into the feed solution, the number of big-sized defects on the surface of fibers decreases, the diameter of fibers is reduced, thinner fibers are packed in a more uniform mode, and the distance between the neighboring fibers is smaller, thus leading to smaller pores. In other words, when the density of the material is reduced (Table 1), the resultant material is characterized by a higher content of air layer between well-cured fibers without bonds, smudges, and thicker regions.

The fine structure of the fibrous *Hmi*/PHB materials was studied in more detail by SEM observations. Figure 3 shows the corresponding SEM images. As follows from Figure 3, the fibrous materials based on neat PHB and *Hmi*/PHB with a low content of *Hmi* (below 3 wt.%) are seen to be non-uniform; the fibers contain joints, crimps, and big-sized elements as spindles. The average size of structural imperfections ranges between 20 and 30 µm along the longitudinal direction and between 15–25 µm in the perpendicular direction. With increasing the content of *Hmi* in the *Hmi*/PHB system, the diameter of fibers markedly decreases, and the fibers within nonwoven materials are seen to be more uniform. Noteworthy is that, for the *Hmi*/PHB systems containing 5 wt.% of *Hmi*, the individual fibers are nearly free of structural defects and imperfections. Both roughness on the fiber surface almost completely, and spindle-shaped thickenings disappeared at 5 wt.% content of Hemin.

To gain a deeper insight into the fine structure of the *Hmi*/PHB materials, elemental analysis, and distribution of *Hmi* additive within the fibers, the *Hmi*/PHB materials were studied by the EDX method (Figure 4).

**Figure 3.** SEM micrographs of PHB and *Hmi*/PHB nonwoven materials with different content of *Hmi*: (**A**) neat PHB, (**B**) *Hmi*/PHB with 1 wt.% of *Hmi*, (**C**) *Hmi*/PHB with 3 wt.%, (**D**) *Hmi*/PHB with 5 wt.% of *Hmi*. **Figure 3.** SEM micrographs of PHB and *Hmi*/PHB nonwoven materials with different content of *Hmi*: (**A**) neat PHB, (**B**) *Hmi*/PHB with 1 wt.% of *Hmi*, (**C**) *Hmi*/PHB with 3 wt.%, (**D**) *Hmi*/PHB with 5 wt.% of *Hmi*.

To gain a deeper insight into the fine structure of the *Hmi*/PHB materials, elemental analysis, and distribution of *Hmi* additive within the fibers, the *Hmi*/PHB materials were studied by the EDX method (Figure 4). According to the EDX analysis, hemin is incorporated into fibers and uniformly distributed within the polymer structure on the surface and in the bulk. Iron and chlorine atoms serve as probes for the characterization of the distribution of *Hmi* within the fibers as the above atoms are the coordination centers of the tetrapyrrole ring. Noteworthy is that the distribution of iron and chlorine atoms appears to be identical (hence, the EDX micrographs illustrating the presence of iron atoms is presented). In Figure 4, iron atoms are seen as orange dots. In the case of the neat PHB, neither iron nor chlorine atoms are seen.

As was mentioned above, PHB is a semicrystalline polymer, which is composed of amorphous and crystalline phases. The crystalline phase is presented by lamellae with folded chains [52]. Under certain conditions formation of spherulites is also allowed [53]. The amorphous phase exists in the rubbery state (glass transition temperature of PHB is 7 ◦C) The influence of technological methods was very great in the formation of the crystal structure of PHB, since it introduced boundary conditions in the process of formation of the

supramolecular structure. Thus, the ES method promoted the formation of pass-through oriented macromolecules [54]. *Polymers* **2022**, *14*, x FOR PEER REVIEW 8 of 19

**Figure 4.** EDX elemental analysis of iron atoms in the *Hmi*/PHB samples with different content of hemin: (**A**) 0 wt.%, (**B**) 1 wt.%, (**C**) 3 wt.%, (**D**) 5 wt.%. **Figure 4.** EDX elemental analysis of iron atoms in the *Hmi*/PHB samples with different content of hemin: (**A**) 0 wt.%, (**B**) 1 wt.%, (**C**) 3 wt.%, (**D**) 5 wt.%.

According to the EDX analysis, hemin is incorporated into fibers and uniformly distributed within the polymer structure on the surface and in the bulk. Iron and chlorine atoms serve as probes for the characterization of the distribution of *Hmi* within the fibers as the above atoms are the coordination centers of the tetrapyrrole ring. Noteworthy is that the distribution of iron and chlorine atoms appears to be identical (hence, the EDX micrographs illustrating the presence of iron atoms is presented). In Figure 4, iron atoms According to the DSC tests, the degree of crystallinity of PHB is ~60–65%. Hence, the content of the amorphous phase in the rubbery state is ~35–40% [55]. As follows from Figure 5, as the content of *Hmi* in the *Hmi*/PHB system increases to 5 wt.%, the degree of crystallinity decreases by 12%. Hence, the presence of *Hmi* has a certain effect on the crystallization of PHB from the solution during the ES process.

are seen as orange dots. In the case of the neat PHB, neither iron nor chlorine atoms are

As was mentioned above, PHB is a semicrystalline polymer, which is composed of amorphous and crystalline phases. The crystalline phase is presented by lamellae with

seen.

tallization of PHB from the solution during the ES process.

through oriented macromolecules [54].

**Figure 5.** Degree of crystallinity versus the content of *Hmi* in the *Hmi*/PHB system. **Figure 5.** Degree of crystallinity versus the content of *Hmi*in the *Hmi*/PHB system.

Of special interest are the profiles of the melting DSC runs curve for the *Hmi*/PHB samples with the different content of *Hmi*. As the content of *Hmi* in the *Hmi*/PHB materials increases, the melting peak upon the first heating run (Figure 6A), becomes more symmetric, the melting temperature remains virtually unchanged, and the degree of crystallinity decreases. Hence, this evidence makes it possible to conclude that, in the presence of *Hmi*, crystallization of PHB proceeds at a slower rate and in a more uniform manner thus leading to the formation of a more regular and well-organized structure. Noteworthy is that, for the samples with a higher *Hmi* content, the low-temperature shoulder in the temperature interval of 150–160 °C becomes less pronounced, thus indicating the lower content of poorly organized crystallites and defects. On the second heating run, the DSC scans show double melting peaks, and this behavior is typical of PHB and corresponds to the melting-recrystallization-remelting process. The melting temperature of initial crystals is observed at 160 °C, whereas the melting temperature after recrystallization is about 174 °C. Of special interest are the profiles of the melting DSC runs curve for the *Hmi*/PHB samples with the different content of *Hmi*. As the content of *Hmi* in the *Hmi*/PHB materials increases, the melting peak upon the first heating run (Figure 6A), becomes more symmetric, the melting temperature remains virtually unchanged, and the degree of crystallinity decreases. Hence, this evidence makes it possible to conclude that, in the presence of *Hmi*, crystallization of PHB proceeds at a slower rate and in a more uniform manner thus leading to the formation of a more regular and well-organized structure. Noteworthy is that, for the samples with a higher *Hmi* content, the low-temperature shoulder in the temperature interval of 150–160 ◦C becomes less pronounced, thus indicating the lower content of poorly organized crystallites and defects. On the second heating run, the DSC scans show double melting peaks, and this behavior is typical of PHB and corresponds to the melting-recrystallization-remelting process. The melting temperature of initial crystals is observed at 160 ◦C, whereas the melting temperature after recrystallization is about 174 ◦C. *Polymers* **2022**, *14*, x FOR PEER REVIEW 10 of 19

folded chains [52]. Under certain conditions formation of spherulites is also allowed [53]. The amorphous phase exists in the rubbery state (glass transition temperature of PHB is 7 °C) The influence of technological methods was very great in the formation of the crystal structure of PHB, since it introduced boundary conditions in the process of formation of the supramolecular structure. Thus, the ES method promoted the formation of pass-

According to the DSC tests, the degree of crystallinity of PHB is ~60–65%. Hence, the content of the amorphous phase in the rubbery state is ~35–40% [55]. As follows from Figure 5, as the content of *Hmi* in the *Hmi*/PHB system increases to 5 wt.%, the degree of crystallinity decreases by 12%. Hence, the presence of *Hmi* has a certain effect on the crys-

**Figure 6.** The DSC scans of Hmi/PHB with different content of Hmi: first heating run (**A**), second heating run (**B**). **Figure 6.** The DSC scans of Hmi/PHB with different content of Hmi: first heating run (**A**), second heating run (**B**).

The Hmi/PHB materials were studied by the EPR method to gain a deeper insight into the inner structure of the amorphous regions of PHB. Noteworthy is that the spin

PHB whereas the crystalline regions contain no TEMPO. Hence, the analysis of the corresponding EPR spectra provides information concerning the mobility of spin radicals within the amorphous regions. Figure 7 shows the correlation times and concentration of TEMPO radicals plotted against the content of *Hmi* in the samples. As follows from Figure 7A, as the content of *Hmi* in the *Hmi*/PHB samples increases, the correlation time tends to decrease, and this decrease suggests the mobility of the TEMPO spin probe is reduced. The content of the spin probe (Figure 7B) appears to be independent of the *Hmi* content in the *Hmi*/PHB samples [54]. In other words, even when the concentration of the probe rad-

**Figure 7.** Correlation times (**A**) and concentration of TEMPO (**B**) versus the content of *Hmi* in the

Marked changes in supramolecular structure and morphology of the *Hmi*/PHB lead to changes in the physical-mechanical properties of the resultant materials (Table 2). For the *Hmi*/PHB containing 5 wt.% of *Hmi*, tensile strength increases up to 5.5 MPa, and this

ical is the same, the mobility of the spin probes decreases.

*Hmi*/PHB material.

heating run (**B**).

The Hmi/PHB materials were studied by the EPR method to gain a deeper insight into the inner structure of the amorphous regions of PHB. Noteworthy is that the spin probe radical TEMPO is localized only within the amorphous phase of semicrystalline PHB whereas the crystalline regions contain no TEMPO. Hence, the analysis of the corresponding EPR spectra provides information concerning the mobility of spin radicals within the amorphous regions. Figure 7 shows the correlation times and concentration of TEMPO radicals plotted against the content of *Hmi* in the samples. As follows from Figure 7A, as the content of *Hmi* in the *Hmi*/PHB samples increases, the correlation time tends to decrease, and this decrease suggests the mobility of the TEMPO spin probe is reduced. The content of the spin probe (Figure 7B) appears to be independent of the *Hmi* content in the *Hmi*/PHB samples [54]. In other words, even when the concentration of the probe radical is the same, the mobility of the spin probes decreases. The Hmi/PHB materials were studied by the EPR method to gain a deeper insight into the inner structure of the amorphous regions of PHB. Noteworthy is that the spin probe radical TEMPO is localized only within the amorphous phase of semicrystalline PHB whereas the crystalline regions contain no TEMPO. Hence, the analysis of the corresponding EPR spectra provides information concerning the mobility of spin radicals within the amorphous regions. Figure 7 shows the correlation times and concentration of TEMPO radicals plotted against the content of *Hmi* in the samples. As follows from Figure 7A, as the content of *Hmi* in the *Hmi*/PHB samples increases, the correlation time tends to decrease, and this decrease suggests the mobility of the TEMPO spin probe is reduced. The content of the spin probe (Figure 7B) appears to be independent of the *Hmi* content in the *Hmi*/PHB samples [54]. In other words, even when the concentration of the probe radical is the same, the mobility of the spin probes decreases.

**Figure 6.** The DSC scans of Hmi/PHB with different content of Hmi: first heating run (**A**), second

*Polymers* **2022**, *14*, x FOR PEER REVIEW 10 of 19

**Figure 7.** Correlation times (**A**) and concentration of TEMPO (**B**) versus the content of *Hmi* in the *Hmi*/PHB material. **Figure 7.** Correlation times (**A**) and concentration of TEMPO (**B**) versus the content of *Hmi* in the *Hmi*/PHB material.

Marked changes in supramolecular structure and morphology of the *Hmi*/PHB lead to changes in the physical-mechanical properties of the resultant materials (Table 2). For the *Hmi*/PHB containing 5 wt.% of *Hmi*, tensile strength increases up to 5.5 MPa, and this Marked changes in supramolecular structure and morphology of the *Hmi*/PHB lead to changes in the physical-mechanical properties of the resultant materials (Table 2). For the *Hmi*/PHB containing 5 wt.% of *Hmi*, tensile strength increases up to 5.5 MPa, and this value is three times higher than that of the PHB nonwoven materials. Elongation at break is equal to 6% (1.7 times higher than that of the initial PHB).

**Table 2.** Mechanical characteristics of the *Hmi*/PHB samples with different content of *Hmi*.

## *3.2. Outdoor Exposure of Hemin/PHB Fibrous Materials*

The stability of the biodegradable *Hmi*/PHB materials containing 5 wt.% of *Hmi* was studied under outdoor exposure. The microscopic images of the test samples before and after weathering experiments are shown in Figure 8. In contrast to the PHB nonwovens

*Polymers* **2022**, *14*, x FOR PEER REVIEW 11 of 19

is equal to 6% (1.7 times higher than that of the initial PHB).

*3.2. Outdoor exposure of Hemin/PHB Fibrous Materials*

**Content of** *Hmi* **%**

containing porphyrins with metal ions when degradation within 60–90 days is complete, the *Hmi*/PHB materials remain nearly intact: the material retained its integrity and shape. containing porphyrins with metal ions when degradation within 60–90 days is complete, the *Hmi*/PHB materials remain nearly intact: the material retained its integrity and shape.

The stability of the biodegradable *Hmi*/PHB materials containing 5 wt.% of *Hmi* was studied under outdoor exposure. The microscopic images of the test samples before and after weathering experiments are shown in Figure 8. In contrast to the PHB nonwovens

value is three times higher than that of the PHB nonwoven materials. Elongation at break

**Elongation at Break (Mean ± SD = 0.2, n = 10) %**

**Table 2.** Mechanical characteristics of the *Hmi*/PHB samples with different content of *Hmi*.

**Tensile Strength (Mean ± SD = 0.05, n = 10) MPa**

0 1.7 3.6 1 0.7 4.7 3 1.9 4.7 5 5.5 6.1

**Figure 8.** Snapshots and SEM images of the *Hmi*/PHB samples before (**A**) and after (**B**) outdoor exposure. **Figure 8.** Snapshots and SEM images of the *Hmi*/PHB samples before (**A**) and after (**B**) outdoor exposure.

The corresponding FTIR spectra of the initial and weathered samples are seen to be nearly the same and reveal no presence of the accumulated products of the degradation process (Figure 9). The corresponding FTIR spectra of the initial and weathered samples are seen to be nearly the same and reveal no presence of the accumulated products of the degradation process (Figure 9). *Polymers* **2022**, *14*, x FOR PEER REVIEW 12 of 19

**Figure 9.** The FTIR spectra of pristine (blue line) and weathered (red line) *Hmi*/PHB samples. **Figure 9.** The FTIR spectra of pristine (blue line) and weathered (red line) *Hmi*/PHB samples.

The FTIR spectrum shows a slight decrease in the intensity of the peaks at 1721 cm−1 (C=O group), 1278 cm−1 (CH3 group), and 1052 cm−1, (C-O-C group). In the region of 900– 700 and 1700–1500 cm−1, noticeable changes are observed. This fact suggests the breakage

Even though the weathered *Hmi*/PHB samples experience minor biodegradation, their supramolecular structure after weathering is changed. The results of the DSC analysis are presented in Figure 10. The shape of the melting peak is markedly changed, thus

**Figure 10.** DSC scans of hemin/PHB (5 wt.%) nanocomposites before (yellow line) and after weath-

ering (blue line): (**A**) first heating run and (**B**) second heating run.

indicating the occurrence of certain processes during the exposure period.

The FTIR spectrum shows a slight decrease in the intensity of the peaks at 1721 cm−1 (C=O group), 1278 cm−1 (CH3 group), and 1052 cm−1 , (C-O-C group). In the region of 900–700 and 1700–1500 cm−1 , noticeable changes are observed. This fact suggests the breakage of C-C bonds and the accumulation of OH groups. The FTIR spectrum shows a slight decrease in the intensity of the peaks at 1721 cm−1 (C=O group), 1278 cm−1 (CH3 group), and 1052 cm−1, (C-O-C group). In the region of 900– 700 and 1700–1500 cm−1, noticeable changes are observed. This fact suggests the breakage of C-C bonds and the accumulation of OH groups. Even though the weathered *Hmi*/PHB samples experience minor biodegradation,

**Figure 9.** The FTIR spectra of pristine (blue line) and weathered (red line) *Hmi*/PHB samples.

Even though the weathered *Hmi*/PHB samples experience minor biodegradation, their supramolecular structure after weathering is changed. The results of the DSC analysis are presented in Figure 10. The shape of the melting peak is markedly changed, thus indicating the occurrence of certain processes during the exposure period. their supramolecular structure after weathering is changed. The results of the DSC analysis are presented in Figure 10. The shape of the melting peak is markedly changed, thus indicating the occurrence of certain processes during the exposure period.

*Polymers* **2022**, *14*, x FOR PEER REVIEW 12 of 19

**Figure 10.** DSC scans of hemin/PHB (5 wt.%) nanocomposites before (yellow line) and after weathering (blue line): (**A**) first heating run and (**B**) second heating run. **Figure 10.** DSC scans of hemin/PHB (5 wt.%) nanocomposites before (yellow line) and after weathering (blue line): (**A**) first heating run and (**B**) second heating run.

The heat of fusion (∆*H*) and melting temperature (Tm) of intact and weathered PHB and *Hmi*/PHB (5 wt.%) nanocomposites are listed in Table 3.

**Table 3.** Thermophysical characteristics (Tm—the melting temperature, ∆*H*—the heat of fusion) of weathered PHB and *Hmi*/PHB nanocomposites (90 days of exposure).

The corresponding stress-strain curves of the *Hmi*/PHB containing 5 wt.% of *Hmi* before and after weathering are shown in Figure 11.

As follows from Figure 11, upon weathering, the mechanical characteristics (tensile stress and elongation at break) are markedly deteriorated. Information on tensile strength and elongation at break is presented in Table 4.

Upon weathering, the tensile strength of the *Hmi*/PHB materials decreases by a factor of 2, and elongation at break decreases from 6.1 down to 4.2%. However, even after outdoor exposure, the *Hmi*/PHB materials (in contrast to the weathering of the neat PHB) preserve their stability and high mechanical properties. Reduced mechanical properties suggest the occurrence of certain irreversible destruction processes promoted by the photosensitizing activity of the tetrapyrrole complex of Hemin. However, the samples preserve good mechanical properties, and both strength and elongation are higher than those of the initial PHB. Slightly increased strength characteristics of PHB are known to be provided by oxidation and UV irradiation.

and *Hmi*/PHB (5 wt.%) nanocomposites are listed in Table 3.

weathered PHB and *Hmi*/PHB nanocomposites (90 days of exposure).

**Days**

**Sample Exposure Time**

fore and after weathering are shown in Figure 11.

**Figure 11.** Strain-stress curves of the initial (blue line) and weathered (yellow line) samples. **Figure 11.** Strain-stress curves of the initial (blue line) and weathered (yellow line) samples.

As follows from Figure 11, upon weathering, the mechanical characteristics (tensile stress and elongation at break) are markedly deteriorated. Information on tensile strength **Table 4.** Mechanical characteristics of the neat and weathered *Hmi*/PHB materials (90 days of exposure).

The heat of fusion (ΔH) and melting temperature (Tm) of intact and weathered PHB

**First Heating Run Second Heating Run Tm, °C ΔH, J/g Tm, °C ΔH, J/g**

**Table 3.** Thermophysical characteristics (Tm—the melting temperature, ΔH—the heat of fusion) of

*Hmi*/PHB (5 wt.%) 0 174.5 88.7 171.0 87.9 *Hmi*/PHB (5 wt.%) 90 171.3 82.8 166.9 80.2

PHB 0 174.1 73.6 170.4 79.4 PHB 90 170.8 63.2 167.2 64.9

The corresponding stress-strain curves of the *Hmi*/PHB containing 5 wt.% of *Hmi* be-

5 90 1.9 4.5

#### 0 0 1.7 3.6 **4. Discussion**

0 90 2.7 4.2 Upon weathering, the tensile strength of the *Hmi*/PHB materials decreases by a factor of 2, and elongation at break decreases from 6.1 down to 4.2%. However, even after Nonwoven materials are characterized by three hierarchical levels of structural organization: macroscopic (the system in whole), mesoscopic (fiber contact area), and microscopic (structure of individual fibers) levels [56], and structural organization at all levels is controlled by the regime of electrospinning and composition of feed solution. For the *Hmi*/PHB composite, the introduction of hemin molecules to the feed mixture appears to exert a marked effect on the structural organization of *Hmi*/PHB at all levels.

At the stage of the electrospinning, the introduction of an even concentration of *Hmi* (1, 3, 5 wt.%) to the *Hmi*/PHB cocktail provides a sustainable regime of ES and leads to the formation of regular nonwoven fibrous material with a uniform distribution of *Hmi* within the amorphous phase between crystallites [51] as proved by the EDX method. The role of *Hmi* is likely to be associated with an increased electrical conductivity of the feed solution due to the presence of the metal ion as the coordination center in the *Hmi* molecule and the viscosity of the solution [57]. As a result, the electrospinning conditions are improved, and the performance of ES is enhanced. As the content of *Hmi* in the *Hmi*/PHB solution increases (from 1 to 5 wt.%), the fibrous material becomes more uniform, and the thickness of regular fibers without defects (smudges, spindles, etc.) decreases down to 1.77 µm. Noteworthy is that the fibers are well developed upon the Taylor cone formation stage [50]. In this case, due to the well-organized movement of a jet upon ES, the porosity of the *Hmi*/PHB material increases from 80 (pristine PHB) to 94% (*Hmi*/PHB containing 5 wt.%). In other words, at the high content of *Hmi* (Table 1), the *Hmi*/PHB material is characterized by a high percentage of the air layer between well-cured fibers.

Concerning the effect of hemin on the morphology of the *Hmi*/PHB materials, the SEM observations (Figure 3) prove that the quality and uniform character of the resultant *Hmi*/PHB nonwoven materials are improved due to the presence of big-sized *Hmi* molecules in the intermolecular space of PHB; as a result, the content of amorphous regions increases, thus increasing the rate of solvent desorption. Evidently, the above changes are provided by the supramolecular structure of PHB [44]. Electrospinning is known to be accompanied by a high deformation rate of the feed polymer solution (103 s−1 ) [51], thus governing the orientation of polymer chains [58]. Upon electrospinning, PHB macromolecules are organized as fibers containing alternating crystalline and amorphous regions [59]. Without the action of external forces, spherulites are formed.

As follows from the DSC scans (Figure 6b), with increasing the content of *Hmi* in the *Hmi*/PHB materials (from 0 to 5 wt.%), the melting peak becomes sharper and the onset of melting is shifted to higher temperatures. This evidence proves the formation of more regular crystallites with regular average dimensions and the absence of poorly organized small crystallites. However, the overall degree of crystallinity decreases (from 65% for neat PHB down to 53% for *Hmi*/PHB containing 5 wt.% of *Hmi*). Hemin is likely to serve as nucleation sites for the formation of uniform PHB crystallites; after crystallization, *Hmi* is accommodated within the amorphous phase of the *Hmi*/PHB materials. According to the EPR observations, as the content of *Hmi* in the *Hmi*/PHB samples increases, correlation time markedly decreases (Figure 8).

The above changes in the morphology of the *Hmi*/PHB materials are associated with the mechanical properties of the resultant materials. The addition of *Hmi* to the feed mixtures provides the formation of regular and uniform fibers and better-organized and well-cured nonwoven materials. The fibers contain no defects as bonds, spindles, smudges, and thinned regions, and their thickness is uniform along their length. Hence, the mechanical characteristics of the *Hmi*/PHB materials are improved.

The stability of the resultant *Hmi*/PHB materials was tested under the conditions of natural weathering (outdoor exposure). In general, weathering of PHB is accompanied by diverse processes, including oxidation, hydrolysis, photooxidation, stress-induced changes, etc., which lead to degradation and cleavage of the polymer chain to shorter fragments [60]. Degradation of biodegradable polymers strongly depends on many factors, including microbial activity, polymer composition, molecular weight, degree of crystallinity, temperature, humidity, pH, UV action, accessibility to oxygen and other oxidizing agents, etc. [61]. According to general knowledge, biodegradation is primarily provided by the biological activity of living organisms and leads to the decomposition of organic compounds into nontoxic products with lower molecular weights. Biodegradation of PHB under anaerobic conditions is accompanied by the formation of carbon dioxide (CO2), water (H2O), and methane; under aerobic conditions, CO2 and H2O are formed [62–64]. Several bacteria and fungi, e.g., *Pseudomonas*, *Actinomadura*, *Penicillium Aspergillus* spp., *Microbispora*, *Saccharomonospora*, *Streptomyces*, *Thermoactinomyces,* and *Bacillusspp.*, can also degrade PHAs both aerobically and anaerobically. Degradation of ultrathin PHB fibers is known to proceed with a high rate as compared with that of PHB films [65] due to the three-dimensional structure and high surface area of nanofibers [65,66].

The results on the outdoor weathering of the *Hmi*/PHB nonwoven materials containing thin fibers seem to be of vital importance for the characterization of their stability and potential use in biomedical applications. As follows from Figure 9, even after outdoor exposure for 90 days, the weathered *Hmi*/PHB samples preserve their shape and integrity. Further tests reveal insignificant fiber cleavage and thinning. The FTIR analysis reveals that degradation proceeds without any accumulation of functional groups or toxic products. As follows from Table 4, this process leads to certain deterioration in the mechanical properties of the *Hmi*/PHB materials: tensile strength decreases by a factor of 2, and elongation at break decreases by a factor of 1.5 (Figure 11). This evidence can be explained by the

degradation-induced changes in the integrity of fibers and supramolecular structure of the *Hmi*/PHB materials due to oxidative and photooxidative processes in the air under the action of the UV light. Earlier, this phenomenon was observed for ozone-induced oxidation of PHB [67].

The DSC data [68] were analyzed for the detailed description of the structural evolution in the crystalline phase as evidenced by marked changes in physicomechanical properties.

Within 90 days of exposure, the supramolecular structure (including the crystalline phase) experiences marked changes in the action of environmental factors (water, wind, UV, temperature changes). The slight decrease in the melting temperature of the weathered samples suggests the formation of poorly organized crystallites and a certain decrease in the molecular weight of the samples [69,70]. However, in general, degradation of the *Hmi*/PHB materials is retarded as compared with the neat PHB. Seemingly, the exposure to solar radiation leads to the initiation of oxidative processes on the surface of individual fibers, leading to minor changes in the crystalline structure as evidenced by a well-pronounced shoulder and double melting peaks which are indicative of the formation of imperfect crystallites. Noteworthy is that the overall degree of crystallinity remains unchanged.

Let us also mention the action of *Hmi* as an efficient antibacterial agent with respect to bacteria and fungi [34,71–73]. In the presence of *Hmi*, intensive generation of hydroxyl radical leads to oxidation of cellular compartments, cell wall rupturing, and membrane leakage, thus providing bacterial apoptosis. As a result, the amorphous phase of the *Hmi*/PHB composite materials, which is the most susceptible to biodegradation, appears to be stabilized. Hence, the resultant *Hmi*/PHB composites with high biocompatibility, good mechanical properties, and retarded biodegradation can be used as effective biomedical materials and scaffolds for regenerative medicine and tissue engineering.

#### **5. Conclusions**

This work addresses the preparation of bioinspired *Hmi*/PHB composite materials with high mechanical properties and retarded biodegradation. The addition of *Hmi* to the feed PHB solution (above 3 wt.%) provides a marked improvement in the performance of the electrospinning process and leads to the formation of well-organized and uniform nonwoven materials with high porosity and good mechanical properties.

The work investigated and described the regularities of changes in the PHB-hemin system during outdoor exposure. The obtained results revealed a deeper understanding of the structural organization of PHB macromolecules and their response to oxidative processes. In the case of the commercial launch of the material, practical recommendations for the use of the PHB-hemin system can be carried out. This material retains its stability within 90 days. The *Hmi*/PHB materials can be recommended for their use as biomedical materials for wound disinfection and accelerated healing, and as scaffolds in regenerative medicine and tissue engineering, and this challenging task requires further detailed studies.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/polym14224878/s1, Figure S1: Photo and schematic view of the single-capillary laboratory unit for the electrospinning process; Figure S2: Histogram the fibers' size distribution of PHB and *Hmi*/PHB nonwoven composites.

**Author Contributions:** Conceptualization, O.V.A.; methodology and experimental design, P.M.T.; validation, O.V.A.; investigation A.K.Z., P.M.T., P.V.I., A.A.O., A.Y.Y. and I.A.V.; DSC tests, P.M.T.; writing—original draft preparation, P.M.T.; writing—review and editing, O.V.A.; supervision, O.V.A.; project administration, A.K.Z.; funding acquisition, P.M.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Russian Science Foundation (Agreement No. 22-73- 00038).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The study was carried out using scientific equipment of the Center of Shared Usage «New Materials and Technologies» of Emanuel Institute of Biochemical Physics and the Common Use Centre of Plekhanov Russian University of Economics.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


**Acknowledgements**

First, I wish to thank Sophie Turcot for her support in organizing and maintaining our laboratory facilities as well as her reliability and diligence in her scientific work. I wish to acknowledge Francis Villeneuve for his methodical approach in the multiple rheological tests conducted. Finally, I wish to extend my most sincere gratitude to Amélie Giroux, Patricia Décarie and Vivianne Chagnon of the *Association professionnelle des nutritionnistes experts en dysphagie* for their patientcentered approach to clinical nutrition.


3.3.3. Recovery of MnCO3 from PLS-2 Solution

Using the same principle, iron precipitation assays were carried out with 1.5 SMR of H2O2 for the oxidation of Fe2+ ions to Fe3+ at pH = 4 and ambient temperature. Within this context, nearly 99.9% of Fe in PLS-2 was precipitated. A concentration of 2754 mg Fe/L was recorded in PLS-2 and a concentration of 2.33 mg Fe/L was measured after Fe removal (Table 5). In this step, the concentration

of Co remained stable, while all Cu was removed and the concentration of Ni decreased from 48.4 to 24.1 mg/L.

To obtain a pure manganese solution, zinc precipitation was achieved by adding Na2S (2 SMR of Na2S at pH 4–5). Under these conditions, up to 100% of zinc was precipitated (Table 5) and the obtained ZnS could be added to the first leaching stage for PLS-1 rich in Zn.

The precipitation of MnCO3 was achieved by the addition of Na2CO3. The quantity of Na2CO3 was calculated for 1 SMR, with pH 4.5–5.0 and at ambient temperature. After rinsing, the collected MnCO3 had a purity of 95%.


*Article* **Assessing the Nationwide COVID-19 Risk in Mexico through the Lens of Comorbidity by an XGBoost-Based Logistic Regression Model**

**Sonia Venancio-Guzmán 1 , Alejandro Ivan Aguirre-Salado 1,\*, Carlos Soubervielle-Montalvo 2 and José del Carmen Jiménez-Hernández 1**

**Abstract:** The outbreak of the new COVID-19 disease is a serious health problem that has affected a large part of the world population, especially older adults and people who suffer from a previous comorbidity. In this work, we proposed a classifier model that allows for deciding whether or not a patient might suffer from the COVID-19 disease, considering spatio-temporal variables, physical characteristics of the patients and the presence of previous diseases. We used XGBoost to maximize the likelihood function of the multivariate logistic regression model. The estimated and observed values of percentage occurrence of cases were very similar, and indicated that the proposed model was suitable to predict new cases (AUC = 0.75). The main results revealed that patients without comorbidities are less likely to be COVID-19 positive, unlike people with diabetes, obesity and pneumonia. The distribution function by age group showed that, during the first and second wave of COVID-19, young people aged ≤ 20 were the least affected by the pandemic, while the most affected were people between 20 and 40 years, followed by adults older than 40 years. In the case of the third and fourth wave, there was an increased risk for young individuals (under 20 years), while older adults over 40 years decreased their chances of infection. Estimates of positive COVID cases with both the XGBoost-LR model and the multivariate logistic regression model were used to create maps to visualize the spatial distribution of positive cases across the country. Spatial analysis was carried out to determine, through the data, the main geographical areas where a greater number of positive cases occurred. The results showed that the areas most affected by COVID-19 were in the central and northern regions of Mexico.

**Keywords:** coronavirus; comorbidity; spatial analysis; logistic regression; ROC curve

### **1. Introduction**

The COVID-19 pandemic has affected many countries worldwide. Among these, the United States, India, Brazil, the United Kingdom and Russia lead the list of countries with the highest number of infected [1]. The outbreak of the SARS-CoV-2 virus was declared a pandemic in March 2020, due to its rapid dissemination and its negative effects across various countries [2]. After 18 months, the United States has reported 82,613,620 confirmed cases and 999,842 deaths; India has reported 43,125,370 confirmed cases and 524,260 deaths; Brazil has reported 30,701,900 confirmed cases and 665,216 deaths, among others. In the case of Mexico (5,752,441 confirmed cases and 324,617 deaths), it is in 15th and 4th place regarding the countries with the larger number of positive cases and deaths, respectively [1].

Coronaviruses are a group of viruses that cause respiratory illnesses [3]. The most common symptoms produced by the virus include fever, dry cough, tiredness and shortness of breath, though, in more severe cases, it causes pneumonia or severe acute respiratory syndrome (SARS) that may lead to death. In other cases, some people infected by the virus

**Citation:** Venancio-Guzmán, S.; Aguirre-Salado, A.I.; Soubervielle-Montalvo, C.; Jiménez-Hernández, J.d.C. Assessing the Nationwide COVID-19 Risk in Mexico through the Lens of Comorbidity by an XGBoost-Based Logistic Regression Model. *Int. J. Environ. Res. Public Health* **2022**, *19*, 11992. https://doi.org/10.3390/ ijerph191911992

Academic Editors: Sachiko Kodera and Essam A. Rashed

Received: 12 July 2022 Accepted: 17 September 2022 Published: 22 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**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/).

do not develop any symptom, but they may still infect the rest of the population. According to the World Health Organization (WHO), the coronavirus can be transmitted from an infected person to others through droplets (also called aerosols) that are expelled when coughing, sneezing or speaking, when shaking a sick person's hand, or by touching objects or surfaces contaminated by the virus [4].

The vulnerability of a person to COVID-19 is determined by various factors; among these, the educational level, age, gender, social conditions, sociodemographic conditions, housing conditions, and even psychological and emotional factors. In their recent study, Ref. [5] analyzed the effects of living conditions on morbidity and mortality in the state of Oaxaca. In this study, it was determined that the health of the population was seriously affected by poor living conditions, and the lack of services such as electricity, gas, water, and health educational services. Therefore, several factors have contributed to the exacerbation of the coronavirus disease. Ref. [6] tested the hypothesis that social vulnerability in Mexico contributes considerably to the probability of hospitalization of people sick from COVID-19. To prove this, a cross-sectional study was carried out with public data from the General Directorate of Epidemiology of the Ministry of Health of Mexico using 5-month period data from individuals with coronavirus. The results show that patients with diabetes have 38.4% probability of being hospitalized, and for those patients whose suffer from some other comorbidity, the risk increases to 42.9% (42.2–43.7%).

Several studies have carried out for the evaluation of comorbidities associated with severe and fatal cases of COVID-19. Ref. [7] reviewed 33 systematic studies and 22 metaanalysis, concluding that, of the total cases, 40% had comorbidities, while fatal cases had about 74% had comorbidities. Hypertension, diabetes, and respiratory diseases were the most frequent comorbidities in severe and fatal cases. Ref. [8] founded that about 66.6% of the deaths caused by from COVID-19 were to men, with a mean age of 69.9 years. In addition, a high incidence of hypertension, diabetes, low platelet levels, and chronic cerebrovascular and cardiovascular diseases were observed among the deceased. Ref. [2] carried out a descriptive analysis of COVID-19 cases in Mexico with data obtained through the official website of the Ministry of Health in Mexico. The study included epidemiological, demographic, and clinical characteristics of the patients who were confirmed positive by a real-time RT-PCR test. The results determined that most of the positive cases of COVID-19 occurred in Mexico City, the average age of the patients confirmed with the disease was 46 years old, most of these cases occurred in people between 30 to 59 years, and 58% of all positive cases occurred in men. Regarding deceased patients, it was determined that they had one or more comorbidities, mainly hypertension, diabetes and obesity. Similarly, Ref. [9] determined that the most influential factor in men and women was obesity, followed by diabetes and hypertension, as well as chronic kidney failure only in the case of women. These findings indicate that these comorbidities were associated with the severity of the disease and predispose to more serious complications of COVID-19.

Multivariate logistic regression has been widely used to identify variables associated with risks of coronavirus disease. Ref. [10] adjusted a multivariable logistic regression models with 450 patients from the Massachusetts General Hospital (MGH). They found that, among patients with diabetes, 42.1% was admitted to the intensive care unit, or ICU, 37.1% required mechanical ventilation, and 15.9% died; while, among patients without diabetes, 29.8% was admitted to ICU, 23.2% required mechanical ventilation, and 7.9% died. In their results, they also found that diabetes and obesity were associated with greater odds of ICU admission and mechanical ventilation. Ref. [11] conducted a study with 220 adult patients with confirmed and suspected COVID-19 using multivariate logistic regression and found that older age was an independent risk factor for mortality. In a study with 648 COVID-19-positive patients with a median age of 34 years, using multivariate logistic regression, Ref. [12] found that independent risk factors for critical outcomes among COVID-19 cases include old age, males, cardiac patients, chronic respiratory diseases, and the presence of two or more comorbidities; all had significant *p*-values < 0.05. In order of predicting COVID-19 severity, Ref. [13] conducted a retrospective study with

287 patients of which 36.6% were classified as severe cases and 63.4% as non-severe cases. They used 23 covariates on blood chemistry and obtained an accuracy of 85.2% and an AUC of 0.928. Ref. [14] analyzed the survival probability of patients with COVID-19 in Mexico considering their characteristics and comorbidities. They carried out logistic regression analyses and fitted nonparametric survival curves using the Kaplan–Meier estimator. A greater risk of death was associated with diabetes, obesity, hypertension, age and gender, and at a lower degree to, obstructive pulmonary disease, kidney disease and immunosuppressive diseases. Ref. [15] implemented two logistic regression models to investigate the rate of hospitalization and mortality against other variables. In the analysis, 10,544 records published by the Epidemiological Surveillance System of Respiratory Viral Diseases of the Ministry of Health of Mexico (Ministry of Health, SSA) were used. The results showed that the majority of the positive cases were men, being 54% times more likely to be hospitalized than women. People older than 50 years and with two or more simultaneous chronic diseases were prone to a higher risk of hospitalization and death.

Currently, there is a wide variety of statistical models and machine learning techniques that allow predictions to be made given a set of characteristics. Ensemble learning uses simple models to form a more precise and efficient algorithm. Boosting is an ensemble technique where simpler algorithms are used in a sequential way to take advantage of each original model in order to improve the precision of the subsequent model (Freund and Shapire, 1996). Some assembly models that are based on the Boosting principle include the AdaBoost [16], CatBoost [17], LightGBM [18] and Extreme Gradient Boosting (XGBoost). XGBoost is an algorithm developed as a research project at the University of Washington by [19] at the SIGKDD conference. This model was built by adding simpler models based on regression trees. The optimization algorithm consisted on the gradient descent on an objective function composed of the sum of individual loss functions of each observation. These loss functions measure the distance between an observation and its prediction based on the sum of its estimate in the previous iteration plus a new function that is added sequentially in each iteration. These characteristics allow the global estimation of the model to be scalable.

The XGBoost algorithm has proven its effectiveness and scalability over other methods in recent studies. For example, Ref. [20] implemented XGBoost to predict PM2.5 concentration per hour. The method was compared to several algorithms including random forest, multiple linear regression, decision tree, and support vector machines. Their results showed that this algorithm superpassed all of these methods. Ref. [21] designed a real driving task to extract data to model driving stress dynamics. They built a driver stress management model based on driving behavior, environment and road familiarity, and performed a cluster analysis with K-means to grouping observations of psychological and driver stress data. The XGBoost was used to monitor stress within each group. Performance comparisons revealed that the XGBoost model significantly transcended support vector machines, random forest and gradient boosting decision tree.

In accordance with the current health emergency around the world, it is necessary to build models that can predict the characteristics of people at higher risk of being infected with COVID-19. Therefore, in this research work, we build an XGBoost model to predict when the patient may or may not present the disease, considering risk factors such as having previous diseases. For this, the proposed methodology to address this problem is initially described. Subsequently, the description of the data is made, as well as the main results obtained from the descriptive statistics and the results predicted by the model. Finally, a brief discussion and conclusion is made about the main results obtained, as well as the scope and limitations of the work carried out.

### **2. Materials and Methods**

The population of Mexico is approximately 126 million people, with a proportion of 51.2% and 48.8% of women and men (Population and Housing Census 2020, [22]), with a median of 29 years old. The country has 32 states, of which the state of Mexico has the largest number of inhabitants, while the state of Colima is the least populated. The data used to carry out this study were collected for the entire country.

### *2.1. Available Data*

The data for this study were obtained from the official website of the General Directorate of Epidemiology of the Ministry of Health of Mexico https://www.gob.mx/salud/ documentos/datos-abiertos-152127 (accessed on 10 May 2022). The records obtained from the database range from February 2020 to the end of April 2022, a total of 15,119,419 registered cases, and 40 explanatory variables. The set of variables used is shown in Table 1. The only available geographic data were the name of the municipality. To obtain the geographic coordinates of the municipalities, a shapefile was used. In order to merge the data, 7,347,705 were eliminated from the database, this due to the presence of missing values on the records and the presence of variables that represented individuals with duplicate characteristics, thus leaving a total of 7,771,714 individuals to be considered for analysis; of these, 3,119,851 (40.14%) was confirmed as positive cases and the rest as negative cases.

**Table 1.** Explanatory variables used in the study.

1 https://www.gob.mx/salud/documentos/datos-abiertos-152127 accessed on 10 May 2022; 2 http://www. conabio.gob.mx/informacion/gis/maps/geo/cabmun2kgw.zip accessed on 10 May 2022.

### *2.2. Logistic Regression*

A generalized linear model is a function that relates a set of response variables with some distribution belonging to the exponential family and a set of independent variables that are obtained through a measurement or observation. This type of model allows the generalization of other existing models by varying the distribution of the response variable or the function that relates *E*(*Yi*) = *µi* with the set of covariates [23]. This last function, commonly called the link function, is a monotonic and differentiable function, which satisfies:

$$\lg(\mu\_i) = X^\top \beta\_i \tag{1}$$

where *X* is an *nxp* matrix of explanatory covariates and *β* ∈ R*p* is a *px*1 vector of coefficients of the model. In the case where the response variable follow a Bernoulli distribution, the parameter *µi* = *pi* represents the probability that a success will occur. The link function commonly used in this case is the *logit* function *g*(*µ*) = *g*(*p*) = log *p* 1−*p* . Consequently, the conditional probability of *Y* given *X* can be expressed as:

$$p(Y = y \mid X) = \frac{e^{X \mid \beta}}{1 + e^{X \top} \beta} = \frac{e^{\beta\_0 + \beta\_1 \mathbf{x}\_1 + \beta\_2 \mathbf{x}\_2 + ... + \beta\_p \mathbf{x}\_p}}{1 + e^{\beta\_0 + \beta\_1 \mathbf{x}\_1 + \beta\_2 \mathbf{x}\_2 + ... + \beta\_p \mathbf{x}\_p}} \tag{2}$$

which satisfies (1). This model is known as the logistic regression model.

### *2.3. XGBoost*

In this research, we used an XGBoost-based logistic regression model for the classification problem, which may be denoted XGBoost-LR, where the explanatory variables are the main risk factors described in [7,8], and the response variable is the presence or absence of COVID-19 disease.

Let be a data set: {(*xi* , *yi*) : *i* = 1, . . . , *n*, *xi* ∈ R*p* , *yi* ∈ R}, where *xi* is the vector of characteristics, and *yi* is the response variable. An assembly model based on regression trees makes use of *K* simpler models, to predict *yi* additively:

$$
\widehat{y}\_i = \sum\_{k=1}^K f\_k(x\_i)\_\prime \qquad f\_k \in \mathbf{F}\_\prime
$$

*y*b*i* is the prediction of *yi* and **F** is the class of functions of all possible regression trees;

$$\mathbf{F} = \left\{ f(\mathbf{x}) = w\_{q(\mathbf{x})} \mid \quad q: \mathbb{R}^p \to I, \quad w \in \mathbb{R}^T \right\},\tag{3}$$

*T* is the number of leaf, *I* = {1, 2, . . . , *T*}, *q*(*x*) represents the index *q*(*x*)-th in vector *w*, *wq*(*x*) represents the *q*(*x*)-th component of *w*, *fk* represents an independent regression tree that corresponds to a tree structure *q* with leaf weights *w*.

Regression trees have a score on each of their leaves, *wi* denotes the score on the *i*-th leaf of a tree. The data used in the training of the model are grouped in the leaf nodes. Thus, in the prediction of an example, the decision rules generated by these trees are considered to classify it in the corresponding leaves and calculate the final prediction adding each of the scores in the leaves given by the vector *w* of each tree [19].

The regularized function to minimize, and to learn each *fi* is:

$$\mathcal{L} = \sum\_{i=1}^{n} l(y\_i, \hat{y}\_i) + \sum\_{k} \Omega(f\_k)\_{\prime} \tag{4}$$

*l* represents the loss function between the observed value *yi* and the predicted value *y*b*i* . In addition,

$$
\Omega(f) = \gamma T + \frac{1}{2}\lambda \|w\|^2.
$$

Ω is the regularized term that penalizes the complexity of the model to avoid overfitting of the data, *γ* penalizes the number of leaves or equivalently the complexity of the tree, *T* is the number of leaves in the tree, *λ* is the regularization parameter, and *w* is the score vector in the leaves.

The loss function is considered as the mean square error if the problem is regression, whereas in a binary classification problem, the loss function is the likelihood function obtained from Equation (2):

$$\text{Let } \hat{y}\_i^{(t)} = \bar{\sum}\_{k=1}^{t} f\_k(\mathbf{x}\_i) = \left\{ \sum\_{k=1}^{t-1} f\_k(\mathbf{x}\_i) \right\} + f\_t(\mathbf{x}\_i), \text{ and Equation (4) can be written as: } \mathbf{x}\_i = \mathbf{x}\_i$$

$$\mathcal{L}^{(t)} = \sum\_{i=1}^{n} l(y\_i, \hat{y}\_i^{(t-1)} + f\_t(x\_i)) + \Omega(f\_t). \tag{5}$$

Using the second order Taylor expansion of the loss function and denoting by *Ij* = {*i* | *q*(*xi*) = *j*} for some fixed structure *q*(*x*) to the index set whose examples correspond to that leaf, we can rewrite (5) as:

$$\mathcal{L}^{(t)}\; \;=\;\sum\_{j=1}^{T}\left[ (\sum\_{i\in I\_{j}}g\_{i})w\_{j} + \frac{1}{2}(\sum\_{i\in I\_{j}}h\_{i} + \lambda)w\_{j}^{2}\right] + \gamma T. \tag{6}$$

The internal sums of the Equation (6) over the index set *Ij* of the leaf *j* is due to the fact that the examples present in that leaf have the same score. Subsequently, to find the optimal value of the parameters *wj* in each of the leaves, fixing *q*(*x*) and deriving Equation (6), it was determined that such a value is:

$$w\_j^\* = -\frac{\sum\_{i \in I\_j} g\_i}{\sum\_{i \in I\_j} h\_i + \lambda}. \tag{7}$$

Substituting (7) in (6), we found that the minimum value is reached in:

$$\mathcal{L}^{(t)}(q) = -\frac{1}{2} \sum\_{j=1}^{T} \frac{(\sum\_{i \in I\_j} g\_i)^2}{\sum\_{i \in I\_j} h\_i + \lambda} + \gamma T.$$

In practice, the goal is to optimize one level of the tree at a time, so that branches are added iteratively to a tree. Thus, if a leaf is divided into two leaves and the index set of the examples that are on the left leaf and on the right leaf are *II* and *ID*, respectively, then the reduction in loss after division is given by:

$$L\_{\text{split}} = \frac{1}{2} \left[ \frac{(\sum\_{i \in I\_I} g\_i)^2}{\sum\_{i \in I\_I} h\_i + \lambda} + \frac{(\sum\_{i \in I\_D} g\_i)^2}{\sum\_{i \in I\_D} h\_i + \lambda} - \frac{(\sum\_{i \in I\_I \cup I\_D} g\_i)^2}{\sum\_{i \in I\_I \cup I\_D} h\_i + \lambda} \right] - \gamma. \tag{8}$$

This last expression consists of the scores of the original leaf, the left leaf, the right leaf, and the regularization term. In addition, this formula is the gain in the loss reduction equation; therefore, it is used to evaluate the division with the candidates.

### **3. Results**

### *3.1. Spatiotemporal Distribution of Daily Cases*

The temporal distribution function of daily cases of COVID-19 in Mexico from February 2020 to the end of April 2022 is described in Figure 1. When analyzing the Figure 1, it can be seen that, after the arrival of the SARS-CoV-2 virus in Mexico, the daily cases of infected increased rapidly beginning in the month of May, reaching the peak of the first wave of positive cases in summer 2020 (mid-July and early August). In this first wave, the temporal distribution function considering the gender of the patient indicates that there is a greater presence of masculine cases than feminine. Subsequently, the second wave of positive cases is confirmed during the winter of 2021 (January–February). In this phase, the occurrence of positive cases was higher than in the first wave, the infection had propagated to every state in Mexico. In addition, notice from Figure 2 that the cases in women and men are quite similar, and the trend of male cases is very similar to the trend of female cases. On the other hand, the third and fourth wave of cases occurred during the months of July 2021 and January 2022. Furthermore, the last wave was the largest, with a significantly greater number of cases in men than in women. The distribution function by age group shows that, during the first and second wave, young people aged ≤ 20 were the least affected by the pandemic, while the most affected were people between 20 and 40 years, followed by adults older than 40 years. In the case of the third and fourth wave, there was an increased risk for young individuals (under 20 years), while older adults over 40 years decreased their chances of infection; see Figure 3.

**Figure 1.** Temporal distribution function in terms of probability for daily COVID-19 cases from February 2022 to April 2021. The dotted line represents the date when a maximum was reached in COVID cases.

**Figure 2.** Temporal distribution function in terms of probability for daily COVID-19 cases by gender.

The results obtained from the statistical analysis by geographic location determined that the states with the highest occurrence of COVID-19 cases were the following: in the central zone of Mexico, the predominant states were Mexico City and the State of Mexico; in the west, Jalisco and Guanajuato; in the north of Mexico, these were Tamaulipas, Nuevo León, Sonora, and, in the south, these were Puebla, Veracruz, and Tabasco, see Figure 4. Most of the states least affected by the disease, with the exception of Chiapas, were the least populated, as is the case of Colima, Nayarit and Campeche; an exception is the state of Chiapas. On the other hand, the municipalities in Mexico with the most cases of COVID-19 were Puebla, Querétaro, León, Monterrey and Guadalajara. In the case of Mexico City, the largest number of cases occurred in the most populated municipalities, such as Iztapalapa, Gustavo A. Madero, Álvaro Obregón, Tlalpan and Benito Juarez, see Figure 5.

**Figure 3.** Temporal distribution function in terms of probability for daily COVID-19 cases by age.

**Figure 5.** Municipalities with most COVID-19 cases from February 2020 to April 2022.

Additionally, the maps in Figure 6a,b show the actual frequency for each municipality of positive cases of COVID-19. Figure 6a shows that many entities in Mexico have a high number of incidents, while Figure 6b shows the center of Mexico, it can be seen that the federal district is one of the most affected entities.

**Figure 6.** Spatial distribution of COVID-19 positive cases: (**a**) national context; (**b**) zoom-in to Central Mexico.

### *3.2. Spatiotemporal Distribution of Individuals with Comorbidities*

In Mexico, a large part of the population suffers from a disease that exposes them to a higher risk of developing COVID-19 or of having serious complications once they are sick. Therefore, it is important to explore the probability of cases against the main risk factors. To do this, four groups that were diagnosed positive for COVID-19 were analyzed: (1) individuals who did not present any risk factor, (2) individuals with diabetes, (3) individuals who suffered from obesity and (4) those who presented pneumonia, Figure 7. For the group of individuals without comorbidity, it was determined that the highest number of incidences occurred in patients between the ages of 25 and 35 years, with an average of 37.75 years. In the case of individuals with diabetes, the mean was presented at

57.5 years, while, for obesity, the cases were positively skewed with a mean of 45.62 years, unlike pneumonia, where the distribution of cases was negatively skewed with a mean of 56.8 years.

**Figure 7.** Probability distribution functions for positive COVID-19 cases by comorbidity against age.

The proportion of patients who presented some risk factor of the total set of records is described in Table 2. Figure 8 shows that the main risk factors present in a large part of the Mexican population are hypertension, pneumonia, obesity, diabetes and smoking. On the other hand, the spatial distribution of each of the comorbidities under study associated with the cases that tested positive for COVID-19 can be seen in Figure 8.

Table 3 shows the patients diagnosed as positive for COVID-19, who had some previous comorbidity. The main risk factors in the Mexican population included hypertension, pneumonia, obesity, and diabetes. In addition, only of the positive cases (3,119,851), 55.59% of the patients suffered at least one of the 10 diseases that were considered in the study. When analyzing the gender, it was found that 50.31% and 49.69% of the total positive cases are men and women, respectively. In the group of men, 56.16% had at least one disease, while in the case of women this was 55%.

**Table 2.** Percentage of patients with any comorbidity, out of a total of 7,771,714 observations.

**Figure 8.** Spatial distribution of comorbidities associated with COVID-19 occurrence from February 2020 to April 2022: (**a**–**j**) all the comorbidities and (**k**) captions for states and countries.

**Table 3.** Percentage of COVID-19 patients with any comorbidity, out of a total of 3,119,851 observations.

### *3.3. Classification Analysis with XGBoost-LR and Classic LR*

This study is focused on determining a model capable of predicting whether or not the patient will suffer from the COVID-19 disease, considering a set of explanatory variables. Additionally, it sought to identify the importance of each of these variables in different groups of patients. The study was performed with the XBGoost-LR model described in Section 2.3 using the likelihood of a multivariate logistic regression model as a loss function. The measure of selection that determined the order of importance of variables in the context

of decision trees was the gain, defined in (8). This measure allowed for deciding which attribute should go to a decision node, selecting the one whose gain is greater.

The XGBoost-LR model was built considering 70% of the data as the training set and 30% as the testing set. To measure model quality among a set of candidate models, evaluation metrics were calculated with only the test set. For the final model, we obtained a score known as sensitivity (66.11%) (i.e., recall rate or fraction of true positives), which is the proportion of positive cases that the model correctly identified. The specificity was 70.11% (true negative rate), which represented the proportion of negative cases that were correctly classified by the algorithm. The probability of how close the result of a measurement is to the true value, that is, the accuracy of the model, was 68.5%. In addition, it was determined that the most influential variables in the model (in terms of decision trees, which generate the greatest gain) for the classification were: date, pneumonia, location with latitude, longitude and age, while diseases such as asthma, hypertension, disease Chronic obstructive pulmonary disease and cardiovascular diseases were the least influential in the model classification. The importance of the variables is shown in Figure 9. This type of importance represents the average gain across all splits where feature was used. Therefore, the table shows that the variables age, latitude, longitude, pneumonia and date are the variables with the greatest importance, while that the variables associated with asthma, hypertension, COPD, cardiovascular disease, another disease, diabetes, chronic, kidney, disease, smoking obesity and gender have lower values of importance. Although the measure of importance obtained by these variables is low, this does not mean that the variables do not contribute to the adjustment; rather, it means that the number of cases in which the variable was decisive to form a new branch is lower.

**Figure 9.** Importance of explanatory variables in XGBoost.

The predictive capacity of the built model was analyzed by comparing it with other classification models. To do this, it was compared with the classic multivariate logistic regression model, which allows for predicting a binary response variable [23]. For training the classical logistic regression model, the same training and testing set used with the XGBoost-LR model was employed. The results obtained for the new model were as follows: the sensitivity was 43.63%, the specificity was 70.49%, and the accuracy was 59.69%. Therefore, it can be showed that the XGBoost-LR model outperformed the classic logistic regression model in almost all the assessment metrics, except for specificity, see Table 4.

Alternatively, the predictive ability of the final model can also be viewed graphically using maps; see Figure 10a,b. The prediction of cases by municipality obtained with the XGBoost-LR model considering the entire Mexican territory can be seen in Figure 10a, while the zoom-in to Central Mexico is found in Figure 10b. Figure 10a shows the level of risk present in each of the municipalities, even more so the areas with the highest probability of risk include the north, coast and center of the country. Note that the map is quite similar to the map in Figure 6, so it can be assumed that the model has good predictive capacity on the data. We also elaborate a map with the predictions using the multivariate logistic regression model to contrast the results of the proposed model. The results are shown in Figure 11, in which we observe a significant difference with respect to the map of the observed values shown in Figure 6.

**Figure 10.** Predicted spatial distribution of COVID-19 occurrence from February 2020 to April 2022: (**a**) average actual occurrence with XGBoost-based logistic regression; (**b**) zoom-in to Central México.

**Figure 11.** Predicted spatial distribution of COVID-19 occurrence from February 2020 to April 2022: (**a**) average actual occurrence with multivariate logistic regression; (**b**) zoom-in to Central México.

The XGBoost-LR model was analytically evaluated by calculating the area under the ROC curve called AUC. The AUC value obtained was 0.75, which indicates that, if an individual is randomly selected, then there is a probability of 0.75 that the diagnosis made by the model is correct. Thus, the constructed model is considered appropriate to predict the risk of acquiring the COVID-19 disease. Figure 12a, shows the ROC curve with the performance of the XGBoost model. Specificity evaluates the proportion of real negatives that are correctly identified, while sensitivity describes the proportion of real positives that were correctly identified. In this way, the ROC curve also provides a perspective on the behavior of the predictive potential of the proposed model. In addition, the optimal cut-off point was also determined, which maximizes the difference between sensitivity and 1-specificity (Figure 12b). Figure 12c shows the ROC curve and AUC obtained for the multivariate logistic regression model. The value of the AUC obtained for this model shows that the proposed model is better than the multivariate logistic regression model. In this model, the cut-off was 0.37, which is slightly lower than the cut-off obtained by the XGBoost model.

**Table 4.** Evaluation metrics of the XGBoost-LR and classic multivariate logistic regression models.

(**a**) ROC curve and AUC of the XGBoost-LR model. (**b**) Global optimal cut-off in 0.4 of the of the XGBoost-

LR model.

(**c**) ROC curve and AUC of the logistic regression model.

(**d**) Global optimal cut-off in 0.37 of the logistic regression model.

**Figure 12.** ROC curve (**a**,**c**), area under the AUC curve (**b**,**d**) and optimal cut-off value of the models: XGBoost-LR and logistic regression, respectively, to estimate the risk of suffering from COVID-19.

### *3.4. Predicting New Cases*

Figure 9 shows that pneumonia is the most important variable of the variables associated with previous comorbidities. In addition, the rest of the variables associated with previous comorbidities also contribute to the construction of the model. The reason is because the decision tree considers spatial and temporal variables, and consequently, in certain geographical places, these variables play an important role for prediction. However, because these places are less numerous within the study region, the importance assigned to the variable is considerably lower. In the case of patients with no history, it was observed that the estimated probability of risk is low, so the aforementioned probability does not exceed the optimal threshold of 0.4, for which the model classifies these observations as negative cases for COVID-19. In the case of patients suffering from obesity and diabetes, it was determined that the probability of risk increases slightly compared to individuals who do not present comorbidity; however, this probability is also below the optimal threshold; thus, they are considered as cases negative for the disease. Finally, in the case of pneumonia, the probability of risk increases beyond the optimal threshold after 10 years.

### **4. Discussion**

Our findings based on descriptive statistics determined that the most prevalent comorbidities in the studied population were obesity, diabetes, hypertension and pneumonia. When analyzed variables of location, the data reflected a notable spatial differentiation in the number of incidents of positive and non-positive COVID-19 cases. That is, urban areas in the central, northern and coastal region of the country, presented a high risk of infection compared to less populated places, such as Chiapas, Campeche and Colima, as can be seen in the maps.

Another result of the descriptive analysis was the ability to determine the time periods with a high rate of contagion in Mexico, which occurred shortly after the holiday seasons of each year. This is reasonable given the enormous mobility on those dates. Additionally, when analyzing the distribution function of positive cases of COVID-19 grouped by comorbidity, it was observed that older adults were more likely at least one of the comorbidities considered in the study. Depending on gender, there was almost the same proportion of positive cases in men and women, which suggests that the risk is independent of gender. On the other hand, after using the XGBoost-LR model, we were able to confirm that, indeed, a large part of the Mexican population was affected by the COVID-19 disease.

Our results also revealed that some vulnerable groups have a greater risk, which may be characterized by their health status and geographic location. This fact is of vital importance not only for the study of the new coronavirus disease, but to show the health status of a large part of the Mexican population. In fact, the Mexican Ministry of Health indicated that, after the United States, Mexico ranks second in terms of obesity rate, a disease that is highly related to diabetes, being one of the main causes of death in the country. On the other hand, the analysis of the geographical locations makes it possible to analyze the existing risk in densely populated places, and encourages the relevant authorities to take better measures, for example, to strengthen security protocols, in touristic areas of the country.

Treviño [24] found the presence of the same comorbidities, finding that men suffer from COVID-19 disease more often than women. In addition, he indicated that admission to the Intensive Care Unit was related to the gender patients, as verified by the Chi-square test. Our result probably differs in that, during the first year of the pandemic, there were more positive cases in men than in women; however, during the second year, this could change. On the other hand Ref. [15] performed a logistic regression analysis and determined that men were 1.54 times more likely to be hospitalized than women. Furthermore, both [25] and [9] found that, in addition to diabetes and obesity, high blood pressure and chronic kidney damage can increase mortality in COVID-19 patients. Regarding the analysis by geographic location of each state of Mexico, Refs. [2,26] analyzed the cases of COVID-19 in Mexico. In both studies, the results agreed with our analysis, since the same spatial pattern

was detected throughout the country: State of Mexico, Nuevo León, among others, were the most affected by the pandemic.

Having a predictive model that allows for identifying individuals with particular characteristics is important for taking preventive measures to avoid further deaths. On the other hand, the constant appearance of new variants of the virus worries stakeholders around the world and makes it difficult to pinpoint an end date for the pandemic. Therefore, the authorities must strengthen surveillance and adopt systematic approaches to provide health indications to the entire population, based on each local context.

### **5. Conclusions**

In the present study, positive and non-positive cases of COVID-19 in Mexico were analyzed at the municipal level using data from the Mexican Ministry of Health. A descriptive analysis was performed to explore the spatial and temporal distribution of positive cases, and subsequently an XGBoost-LR model was trained to estimate the risk of infection using the main comorbidities in patients as covariates. According to the results, we confirm that pneumonia, obesity and diabetes were variables of utmost importance to predict new cases. The estimated and observed values of percentage occurrence of cases were very similar, and indicated that the proposed model was suitable to predict new cases (AUC = 0.75). The optimal global cutoff to identify when a case is positive for COVID-19 was 0.41. Our results revealed that patients without comorbidities are less likely to be COVID-19 positive, unlike people with diabetes, obesity and pneumonia with ages close to 56, 46 and 59 years old, respectively. The current worldwide situation regarding COVID-19 is still serious, especially with the emergence of new SARS-CoV-2 virus variants. Therefore, we recommend reviewing and adapting government policies concerning the establishment of preventative measures to avoid the spread of virus and mitigate the pandemic effects.

**Author Contributions:** Conceptualization, S.V.-G., A.I.A.-S., C.S.-M. and J.d.C.J.-H.; Data curation, S.V.-G., A.I.A.-S., C.S.-M. and J.d.C.J.-H.; Formal analysis, S.V.-G., A.I.A.-S., C.S.-M. and J.d.C.J.-H.; Funding acquisition, S.V.-G., A.I.A.-S., C.S.-M. and J.d.C.J.-H.; Investigation, S.V.-G., A.I.A.-S., C.S.-M. and J.d.C.J.-H.; Methodology, S.V.-G., A.I.A.-S., C.S.-M. and J.d.C.J.-H.; Project administration, S.V.-G., A.I.A.-S., C.S.-M. and J.d.C.J.-H.; Resources, S.V.-G., A.I.A.-S., C.S.-M. and J.d.C.J.-H.; Software, S.V.-G., A.I.A.-S., C.S.-M. and J.d.C.J.-H.; Supervision, S.V.-G., A.I.A.-S., C.S.-M. and J.d.C.J.-H.; Validation, S.V.-G. and A.I.A.-S.; Visualization, S.V.-G. and A.I.A.-S.; Writing—original draft, S.V.-G. and A.I.A.-S.; Writing—review & editing, S.V.-G. and A.I.A.-S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by a CONACYT Grant scholarship: 2021–2023, for the doctoral studies of the first author. This research was funded by CONACYT through the grant "Convocatoria de Ciencia Básica y/o Ciencia de Frontera 2022", project ID 320036.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The analyzed data sets can be consulted at the following links: https: //www.gob.mx/salud/documentos/datos-abiertos-152127 accessed on 10 May 2022; http://www. conabio.gob.mx/informacion/gis/maps/geo/cabmun2kgw.zip accessed on 10 May 2022.

**Acknowledgments:** The authors thank the Secretary of Health of Mexico (https://www.gob.mx/ salud/documentos/datos-abiertos-152127, accessed on 10 May 2022) for providing the data used in this research. Special thanks are also given to two anonymous reviewers who shared insightful observations that deeply improved our work.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Abbreviations**

The following abbreviations are used in this manuscript:

### **References**