**Celebrating Applied Sciences Reaches 20,000 Articles Milestone: Feature Papers in Applied Biosciences and Bioengineering Section**

Editors

**Francesco Cappello Magdalena Gorska-Ponikowska**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Francesco Cappello University of Palermo Italy

Magdalena Gorska-Ponikowska Medical University of Gdansk Poland

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: https://www.mdpi.com/journal/applsci/special issues/Celebrating Biosciences Section).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-4427-4 (Hbk) ISBN 978-3-0365-4428-1 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**

## **About the Editors**

#### **Francesco Cappello**

Francesco Cappello was born in Palermo, Italy, in 1973. In 1997, he completed his degree in Medicine and Surgery with honors and was awarded with Best Thesis prize. In 2002, he became a specialist in Pathological Anatomy with the maximum score. Since 2016, he has served as Full Professor at the University of Palermo, Palermo, Italy, where he teaches Human Anatomy and Histology. His current honorific appointments include the following: Scientific Director of the Euro-Mediterranean Institute of Science and Technology, Palermo, Italy; Associate Member of the Faculty at the Neuroscience Graduate Program, University of Texas Medical Branch, Galveston (TX), USA; Adjunct Associate Professor at the Department of Biology, Temple University, Philadelphia (PA), USA; President of the Italian Society of Experimental Biology; Senior Fellow of the Cell Stress Society International. His research interests include cell differentiation, tissue homeostasis, organ regeneration, cell stress and chaperones, and nanovesicles.

#### **Magdalena Gorska-Ponikowska**

Magdalena Gorska-Ponikowska (born in 1984), Associate Professor, is the Head of the Department of Medical Chemistry, Medical University of Gdansk, Poland. She has been a Visiting Researcher at the University of Stuttgart, Germany, since 2017 and has also served at the Institute 'Istituto Euro Mediterraneo di Scienza' (Palermo, Italy) since her appointment in 2018. She is a creator and co-owner of the cosmetic brand Skin Science, contributing with her expert authorship. Prof. Magdalena Gorska-Ponikowska graduated from the Faculty of Pharmacy, Medical University ´ of Gdansk, Poland. She subsequently obtained her habilitation degree in 2018. She has participated in numerous international scientific research internships, including at University of Palermo (Italy), University of Stuttgart (Germany), and University of Alberta (Edmonton, Canada). Prof. Magdalena Gorska-Ponikowska was awarded a scholarship from the Minister of Science and Higher Education ´ (2017) and granted the title 'Muse of the City of Sopot' (Poland) in the field of regional scientific activity (2017). Her research areas involve medical chemistry, biochemistry, molecular biology, oncology, cosmetology, nanotechnology, and neurobiology. She is an author and co-author of 74 peer-reviewed articles in international journals. ORCID: 0000-0002-7366-8429.

### *Editorial* **Special Issue "Celebrating Applied Sciences Reaches 20,000 Articles Milestone: Feature Papers in Applied Biosciences and Bioengineering Section"**

**Magdalena Gorska-Ponikowska 1,2,3,\* and Francesco Cappello 3,4,\***

This Special Issue celebrates the publication of 20,000 articles in *Applied Sciences*. This Special Issue intended to collect papers featuring important and recent developments or achievements in biosciences and bioengineering, with a special emphasis on recently discovered techniques or applications. Thanks to it, we managed to gather the interdisciplinary papers, including the broad spectrum of the following topics: clinical studies, advanced diagnostics, biochemistry, bioinformatics, biomaterials, biomechanics, biomedicine, biotechnology, clinical engineering, drug delivery, microbiology, life science, neuroscience, oncobiology, physiology, radiology, and tissue engineering.

In vivo clinical studies are the most important determination of drugs' pharmacokinetics and efficacy [1–4]. In an intensive care unit clinical study, the pharmacokinetics of levetiracetam, a second-generation antiepileptic drug has been evaluated [1]. Interestingly, a controlled clinical trial supervised by Lai et al. proved that combining yoga with rehabilitation has the potential to improve depressive disorders [2], while Khan and Reilly proposed suggestions for enhanced confidence for randomized controlled trials in protective treatments against endothelial glycocalyx degradation in surgery [3]. The effects of essential oils and other substances derived from the Lamiaceae family plants as adjuvants for the treatment of periodontitis have also been discussed [4].

A number of works concerned the very important topic of searching for biomarkers of pathogenesis and the progression of pathologies, such as cancer. Heat shock proteins (HSPs) are ubiquitously expressed housekeeping chaperones responsible for maintaining homeostasis of the organisms and can be considered as physiologically expressed biomarkers of cancer, e.g., leukemia, or different pathologies including gastric diseases [5–9]. The biomarkers may be localized extracellular or extracellularly excreted [8,9]. The molecular mechanism of asthma and COPD based on extracellular nanovesicles and their putative use in therapy is discussed by Fucarino et al. [8], while Alberti et al. describe tumor-secreted extracellular vesicles as the main mediators of cell–cell communication, permitting cells to exchange proteins, lipids, and metabolites under varying pathophysiological conditions [9].

When it is not possible to perform in vivo tests, establishing an in vitro-based system that can realistically simulate in vivo conditions is desirable [10–14]. A new cell culture method by combining fluoropolymers and dot-patterned extracellular matrix substrates to achieve spheroids has thus been successfully developed [10]. A new low-cost and simple-to-use method for the determination of free biothiols in biological fluids has been also proposed [11]. Moreover, a bioinformatics-based method, which introduces thermodynamic measures and topological characteristics aimed to identify potential drug targets for pharmacoresistant epileptic patients has been established [13]. Di Bella et al. show a relative accuracy, sensitivity, and specificity of 100% for *Salmonella* spp. detection

**Citation:** Gorska-Ponikowska, M.; Cappello, F. Special Issue "Celebrating Applied Sciences Reaches 20,000 Articles Milestone: Feature Papers in Applied Biosciences and Bioengineering Section". *Appl. Sci.* **2022**, *12*, 3978. https://doi.org/10.3390/ app12083978

Received: 6 April 2022 Accepted: 11 April 2022 Published: 14 April 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/).

and identification in comparison with the reference method ISO 6579-1:20 [14]. Interestingly, it was determined that the temperature of storage up to 7 months does not significantly affect the antioxidant properties of elderberry (*Sambucus nigra* L.) juice, which is highly rich in polyphenols, particularly flavonoids [15]. The nutritional characteristics of *Halimione portulacoides* (L.) has been widely described [16].

**Author Contributions:** Conceptualization M.G.-P. & F.C.; writing—original draft preparation, M.G.-P. & F.C. writing—review and editing, M.G.-P. & F.C.; supervision, M.G.-P. & F.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** M.G.-P. kindly acknowledge ST46 funding from Medical University of Gdansk.

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

#### **References**

### *Article* **The Pharmacokinetics of Levetiracetam in Critically Ill Adult Patients: An Intensive Care Unit Clinical Study**

**Sophia-Liberty Markantonis 1, Nikolaos Markou 2, Apostolos Karagkounis 1, Dionysia Koutrafouri 2, Helen Stefanatou 2, Rania Kousovista 3 and Vangelis Karalis 1,\***

<EMAIL_ADDRESS>(N.M<EMAIL_ADDRESS>(D.K<EMAIL_ADDRESS>(H.S.) 3 Department of Mathematics, University of Crete, Heraklion, 70013 Crete, Greece<EMAIL_ADDRESS>
**\*** Correspondence<EMAIL_ADDRESS>Tel.: +30-210-7274267

**Abstract:** The aim of this study was to investigate levetiracetam pharmacokinetics in critically ill adult intensive care patients and to identify pathophysiological factors affecting its kinetics. Fourteen critically ill patients in an intensive care unit were enrolled in the study and received intravenous levetiracetam. Blood samples were collected at specific time points to determine the levetiracetam pharmacokinetics. Patient characteristics such as renal function, demographics, disease severity, organ dysfunction, and biochemical laboratory tests were evaluated for their influence on the kinetics of levetiracetam. Estimated glomerular filtration rate (eGFR) had a statistically significant (*p* = 0.001) effect on levetiracetam clearance. None of the other patient characteristics had a statistically significant effect on the pharmacokinetics. Simulations of dosing regimens revealed that even typically administered doses of levetiracetam may result in significantly increased concentrations and risk of drug toxicity in patients with impaired renal function. The Acute Physiology and Chronic Health Evaluation II (APACHE II) score differed significantly among the three groups with different epileptic activity (*p* = 0.034). The same groups also differed in terms of renal function (*p* = 0.031). Renal dysfunction should be considered when designing levetiracetam dosage. Patients with a low APACHE II score had the lowest risk of experiencing epileptic seizures.

**Keywords:** levetiracetam; critically ill patients; population pharmacokinetic modeling; simulated dosage schemes; intensive care unit

### **1. Introduction**

Levetiracetam, a second-generation antiepileptic drug (AED) with a unique mechanism of action, is approved by both the U.S. Food and Drug Administration and the European Medicines Agency (EMA) as an add-on therapy for the treatment of partial seizures, myoclonic seizures, and primary generalized tonic-clonic seizures. In Europe, it is also the only AED approved for the treatment of partial seizures in adults and adolescents aged 16 years and older [1].

In critical care, levetiracetam has been studied for the treatment of status epilepticus, tumor-related seizures, and seizures following subarachnoid or intracerebral hemorrhage, trauma, and stroke [2–7]. The pharmacokinetic and pharmacodynamic properties of this antiepileptic drug make it an advantageous option in the care of patients in the ICU. Compared to other AEDs, the pharmacokinetic profile of levetiracetam is favorable after both oral and intravenous administration (IV) [8,9]. After oral administration, there is rapid and nearly complete absorption (95%) with dose-proportional pharmacokinetics, low protein binding (10%), and low within-subject variability. It is excreted primarily renally by glomerular filtration with partial tubular reabsorption and has a plasma elimination half-life of approximately 6–8 h in adults. There are no known clinically significant drug

**Citation:** Markantonis, S.-L.; Markou, N.; Karagkounis, A.; Koutrafouri, D.; Stefanatou, H.; Kousovista, R.; Karalis, V. The Pharmacokinetics of Levetiracetam in Critically Ill Adult Patients: An Intensive Care Unit Clinical Study. *Appl. Sci.* **2022**, *12*, 1208. https://doi.org/10.3390/ app12031208

Academic Editors: Francesco Cappello and Magdalena Gorska-Ponikowska

Received: 30 November 2021 Accepted: 23 January 2022 Published: 24 January 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/).

interactions, and it is not associated with the common hemodynamic or cardiovascular side effects of other AEDs [8,10–12].

A clear association between levetiracetam serum levels and therapeutic or toxic effects has not been demonstrated [10,13]. In accordance with the International League Against Epilepsy (ILAE) guidelines for medication adherence monitoring, overdose and dose adjustment, a therapeutic range of 12–46 mg/L has been recommended [10,14,15]. However, different levetiracetam levels outside the reference range in different populations have been reported in the literature. Pregnant women, elderly and pediatric patients, and patients with end-stage renal failure or those requiring renal replacement therapy exhibit altered pharmacokinetics [16–27].

Relatively few studies in critically ill ICU patients have reported increased drug clearance (greater than 120–160 mL/min/1.73 m2) and treatment failure for seizures [21,28–31]. Based on empirical evidence, therapeutic drug monitoring for levetiracetam has been classified as 'potentially useful' but increasing evidence that altered pharmacokinetics of levetiracetam may be observed in critically ill patients has increased the need for therapeutic drug monitoring and pharmacokinetic modeling to determine optimal dosing regimens to maximize therapeutic effect and limit toxicity [22].

The aim of this study was to investigate the pharmacokinetics of levetiracetam in critically ill intensive care patients using nonlinear mixed-effect modeling approaches. The latter allow the identification of important pathophysiological factors (such as renal function) that could influence the kinetics of levetiracetam. Special emphasis was placed on the selection of the most suitable index for renal function and for this purpose several indices have been measured and investigated. Most of the patients participating in the study had from mild loss to severe loss of kidney function which allowed this exploration. Simulations were performed to assess the impact of important patient characteristics on the kinetics of levetiracetam. In addition, the potential influence of disease severity and renal function on epileptic status was investigated.

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

#### *2.1. Patients and Data Collection*

The study was conducted in the intensive care unit of a tertiary hospital (Latsio Burn Center, Thriassio General Hospital of Elefsina, Attica, Greece). Fourteen critically ill patients were enrolled in the study and all of them received intravenous (IV) levetiracetam (Keppra®) either for prophylaxis or treatment of epilepsy. The study protocol was reviewed by the hospital Scientific Committee and written informed consent was obtained from all patients or their legal representatives before enrolment in the study. The study was conducted in accordance with the International Conference on Harmonization guidelines for good clinical practice and the Declaration of Helsinki, while all patient data were collected and processed in full compliance with the EU General Data Protection Regulation 2016/679 and EU Directive 2016/680 of the European Parliament [32–35].

Patient data included information on levetiracetam administration and ICU admission, demographics, biochemical parameters (e.g., liver enzymes, urea, albumin), fluid balance (in mL) and organ function expressed by the Sequential Organ Failure Assessment Score (SOFA) and the Acute Physiology and Chronic Health Evaluation II (APACHE II) score [36]. From patients' weight and height, the body mass index (BMI) and body surface area (BSA) were also calculated [37]. The SOFA score is used to track the status of an ICU patient and determine the extent of their organ function or rate of failure. In this study, several expressions of this score were calculated, such as SOFA hemodynamic, SOFA liver, SOFA kidneys, and overall SOFA. The APACHE II is calculated within 24 h of a patient's admission to the ICU using various measurements such as heart rate, respiratory rate, blood pH, and hematocrit. Taking all these factors into account, an integer value from 0 to 71 is calculated, with higher values indicating a more severe condition and a higher risk of death. Epileptic activity was also monitored and classified as negative, possible, and positive. The "positive" group included patients who had epileptic seizures, while

the "negative" group included those who had no epileptic activity and levetiracetam was administered for prophylactic reasons only. The "possible" category included patients who had a clinical profile of epileptic activity, but this was not verified by electroencephalogram.

In this study, particular attention was paid to the role of renal function, as many of the patients in the ICU had impaired renal function. For this reason, various measures of renal function were calculated, such as creatinine clearance over time, Cockcroft–Gault creatinine clearance, and measures of estimated renal function (eGFR) such as MDRD eGFR (Modification of Diet in Renal Disease study group), CKD-EPI eGFR (Chronic Kidney Disease epidemiology collaboration), Jelliffe eGFR in mL/min/1.73 m2 [38–40]. Patients in this study were divided into two subgroups according to their renal function; patients with stable renal function (group A) and patients with acute renal impairment (group B) (Table 1). Levetiracetam treatment data included dosing regimen, cause of administration, chronic treatment, and loading dose administration.

**Table 1.** Patients' demographic characteristics, levetiracetam dosing regimen, epileptic activity, and subgroup classification.

Key: Epileptic activity: - no activity, ± possible, + epileptic activity; BID, twice daily; eGFR, estimated glomerular filtration rate; Stable renal function, when eGFR > 60 mL/min/1.73 m2.

Levetiracetam was administered diluted in 100 mL N/S 0.9% *w/v* by IV infusion over 15 min, while the loading dose was administered in 250 mL N/S 0.9% *w/v* by IV infusion over 30 min. Blood samples for all participants were taken two days after the start of dosing to ascertain whether the serum levels had achieved a steady state. The specific time points for blood collection were determined according to the pharmacokinetics of levetiracetam to obtain as much information as possible. The schedule for blood collection was 0 h (trough level), 20 min, and 3 h after steady-state administration of levetiracetam. The actual blood sampling times were recorded.

For blood collection and storage of serum or plasma samples the following procedure was applied: collection of 5 mL of whole blood in a coated Vacutainer tube with a red tip (Becton Dickinson). After blood collection, it was allowed to clot undisturbed at room temperature (usually 15–30 min). The clot was removed by centrifugation at 1000–2000 rpm for 10 min. After centrifugation, the serum was immediately divided into two fractions >200 μL and immediately transferred to polypropylene tubes (not gel tubes) where they were stored in a freezer −70 ◦C. For pharmacokinetic analysis, serum concentrations of levetiracetam were determined by high performance liquid chromatography using a UV detector at 205 nm after sample preparation by solid phase extraction using ClinRep® HPLC Complete Kit, levetiracetam (Keppra®) in serum/plasma (order number 15500) (RECIPE Chemicals & Instruments GmbH, Munich, Germany). The lower limit of detection was 0.14 mg/L, the lower limit of quantification was 0.46 mg/L, the upper quantification limit was 104 mg/L, and the recovery was 97–105%. The chromatographic test was performed with the ClinTest® standard solution, calibration with the ClinCal® serum calibrator and quality control with the ClinChek® serum controls level I & II test solutions (order number 15582). The precision of the method was determined by calculating the relative standard deviation at three plasma concentrations during the same analysis (within-day precision) and in triplicate over six analyses (between-day precision). For the entire concentration range, the relative standard errors for both within- and between-assay precision were less than 7.2%.

#### *2.2. Statistical Analysis*

All statistical analyses of dosing schedules and clinical data were performed using IBM SPSS Statistics version 25 (IBM Corporation, Armonk, NY, USA). Statistical comparisons were performed at 5% significance level. For scale variables such as SOFA and eGFR measurements, a normality test was first performed (using the Shapiro–Wilk test) to determine whether the variables followed a normal distribution. In the case of normally distributed data, the independent t-test was performed to compare these variables between the two groups of patients. For comparison of more than two groups, the one-way method ANOVA was used with the post-hoc criteria of least significant difference and Tukey. The non-parametric analogue, Mann–Whitney, was used for variables that deviated from the normal distribution. The relationship between two nominal (or ordinal) measures (e.g., patient group and epilepsy status) was examined using the chi-square test. APACHE II score is actually an ordinal scale, but with many values. In other words, the APACHE II score can be considered as a Likert scale and therefore parametric statistical methods can be applied.

#### *2.3. Population Pharmacokinetic Analysis*

#### 2.3.1. Non-Linear Mixed Effect Modeling

Population pharmacokinetic analysis was implemented in MonolixTM 2020R1 (Lixoft, Orsay, France, Simulation Plus) where individual serum concentration profile data were analyzed using the expectation maximization algorithm of stochastic approximation for nonlinear mixed effects followed by importance test methods. The value of the objective function was calculated using the Monte Carlo method for the final population parameter values.

#### 2.3.2. Structural Models, Error Models, and Covariates

One- and two-compartment models for intravenous infusion with first-order elimination and initial estimates for the parameters were examined [23,25,26,32–34]. A lognormal distribution of pharmacokinetic parameters was assumed, while several residual error models were tested, including constant, proportional, and combined. The parameters were estimated using the stochastic approximation estimation method (SAEM), and the objective function value was determined using the importance sampling Monte Carlo approach at the final population parameter values.

Once the structural model was determined, several covariates were tested. The covariates examined in this study were related to subject-specific characteristics, including sex, age, body weight, height, APACHE II score, timed creatinine clearance, Cockcroft– Gault eGFR, Jelliffe eGFR, CKD-EPI eGFR, MDRD eGFR, serum creatinine, urea, SGOT, and total SOFA. Three expressions of body weight were used: typical body weight, ideal body weight, and adjusted body weight. Linear and lognormal (allometric) models were assessed. Analyses of covariates were performed using stepwise forward selection and backward elimination. Continuous covariates were examined either untransformed or centered around the mean or median value of the covariate. The categorical covariate

examined was sex, while SOFA hemodynamic, SOFA renal, and SOFA liver were ordinal variables. The total SOFA score, which comes from the sum of all individual SOFA scores, can be considered a continuous variable

The Pearson correlation test and one-sided ANOVA were used for continuous and categorical covariates, respectively. The Wald test was used to check whether the covariates could explain the variability of the parameters in the final model. For all pharmacokinetic analyses, the significance threshold was set at 5%.

#### 2.3.3. Model Evaluation

Model selection was based on goodness-of-fit criteria, visual inspection of diagnostic plots, comparison of relative standard deviations of estimated parameters, precision of estimates, and changes in Akaike and Bayesian information criteria and log-likelihood. Visual inspection of goodness-of-fit was performed using plots of observed values versus predicted values for the population and individual weighted residuals versus time. Visual predictive checks (VPCs) were used to assess the predictive performance, stability, and robustness of the model. VPCs were generated from 1000 Monte Carlo simulations and 90% prediction intervals. Normalized prediction distribution errors versus concentration were also used.

#### 2.3.4. Simulations

The final pharmacokinetic model was used to simulate the concentration–time profiles of levetiracetam at different levels of renal function. Three simulated subject groups of 50 subjects each were formed: (a) subjects with normal renal function using the model parameters reported in the literature [23], (b) subjects with 50% clearance limitation, and (c) subjects with the pharmacokinetic characteristics determined in the fourteen patients in this study. The dosing regimen for each subject was identical to that actually used in the ICU (see Table 1). All simulations were performed in Simulx® (MonolixTM 2020R1).

#### **3. Results**

Eleven patients with stable renal function in the ICU had an eGFR estimate greater than 60 mL/min/1.73 m2 and were classified in subgroup "A". The other three patients had acute renal failure and were classified in subgroup "B" (Table 1).

Regarding the reasons for administration of levetiracetam, most patients (71.4%) received levetiracetam for prophylactic purposes and only 28.6% for therapy. Only one patient (i.e., 7.1%) received levetiracetam for chronic treatment. The proportion of patients receiving a loading dose of levetiracetam was 42.9%, while the higher proportion (57.1%) initiated levetiracetam without a loading dose. The reason for levetiracetam administration was either prophylaxis (71.4% of patients) or treatment of seizures (28.6%). The mean eGFR of patients was 94.3 mL/min/1.73 m2 (or equivalent to 5.66 L/h/1.73 m2) and the mean APACHE II score was 17.4.

#### *3.1. Epileptic Activity*

The possible influence of disease severity (expressed by the APACHE II score or the SOFA score) and renal function (expressed by the CKD-EPI eGFR) on epileptic status was investigated. The Shapiro–Wilk test indicated normal distribution of the data and therefore one-way ANOVA was performed to detect a difference in the values of APACHE II and eGFR between the three levels of epileptic activity (positive, negative, possible). It was found that the APACHE II score was significantly different between the three groups (*p* = 0.034). A statistically significant difference was also found for the CKD-EPI eGFR (*p* = 0.031). For the APACHE II score, the absence of epileptic activity (i.e., the "negative" group) differed from the "possible" (*p* = 0.033) and the "positive" epileptic activity (*p* = 0.027). The patients with the lowest APACHE II score belonged to the group without epileptic seizures.

Regarding the CKD-EPI eGFR values, the three epilepsy groups showed significant differences (*p* = 0.031). Similarly, the group of patients without seizures was statistically different from the "possible" (*p* = 0.034) and the "positive" group (*p* = 0.022). Statistical comparisons were made for all other characteristics, but no significant differences were found.

It should be clarified that the abovementioned analysis was not the primary purpose of the study and for this reason no power assessment was made focusing on this. In addition, this assessment only attempted to link disease severity with renal function and not to explore their relationship with levetiracetam pharmacokinetics. In any case, it should be underlined that due to the limited sample size and therefore the number of patients within each group, these findings should be considered as exploratory and further adequately powered studies are necessary to confirm the results.

#### *3.2. Population Pharmacokinetic Model*

Population pharmacokinetic analysis was initially performed separately for each patient group. However, due to the limited number of ICU patients in group B, no model could be derived. Therefore, in order not to waste the information of group B patients, all individuals were pooled into one group to increase the ability to develop a robust pharmacokinetic model, but each patient's individual characteristics (e.g., renal function) were assessed and all available data were used. Furthermore, this pooling was possible because there was no statistically significant difference (*p*-value > 0.05) in levetiracetam concentration values between the two groups.

A one-compartment model with first order elimination was finally selected as the best model for the description of the levetiracetam concentration—time profiles. The estimates of the pharmacokinetic parameters of this best model are given in Table 2. It is of note that the volume of distribution of the central compartment (V) was 47.19 L and the clearance (CL) was 1.45 L/h. Of all the covariates examined for their influence on the kinetics of levetiracetam, only CKD-EPI eGFR was found to be significant for the clearance of levetiracetam. The coefficient of CKD-EPI eGFR with respect to clearance was 0.25, implying that better renal function (in terms of eGFR) leads to higher clearance of levetiracetam. The residual error model that led to the best results was the proportional error model with b = 0.099. The percentage relative standard errors for all estimated parameters had relatively low values. The significant covariate in the final model was the CKD-EPI eGFR (Wald test *p* = 0.00164), which is the Chronic Kidney Disease Epidemiology Collaboration equation, developed in an effort to provide a more precise formula for estimating glomerular filtration rate using serum creatinine, age, sex, and race data.

**Table 2.** Population parameters of the final pharmacokinetic model of levetiracetam.

Key: V, volume of distribution; Cl, clearance; ω\_V, between-subject variability value for V; ω\_Cl; between-subject variability value for Cl; beta\_eGFR, allometric scaling factor for CKD-EPI eGFR (centered around median) on Cl; b, proportional component of the error model.

Figure 1 shows the graphical evaluation of the final model. The visual representation of the predictive power shows that the prediction interval of the developed model includes the experimental concentration data in all cases. Similarly, the good predictive ability of the model is also illustrated in Figure 2, where the population predicted versus observed concentration values are almost linearly correlated (Figure 2a) and the NPDE versus time (Figure 2b) and concentration (Figure 2c) plots show an adequate performance. Additional goodness-of-fit plots are shown in the supplementary material (Figure S1a–c).

**Figure 1.** Visual predictive check plot for the final best levetiracetam pharmacokinetic model. The blue lines refer to the 10th, 50th, 90th percentile of empirical data and the shaded areas refer to the predicted 90% confidence intervals around each zone (10th, 50th, 90th percentiles). A number of 1000 Monte Carlo simulations were used.

The shrinkage of the volume and clearance estimates (parameter values were randomly sampled from the conditional distribution) was adequate and equal to 18.3% and −3.23%, respectively. The relationship between CKD-EPI eGFR and levetiracetam clearance is shown in Figure 3. An almost linear relationship (*levetiracetam clearance = 0.01534*·*(CKD-EPI eGFR)* − *0.21031*) is observed with a correlation coefficient of 0.76. This finding is further evidence that an increase in CKD-EPI eGFR is associated with a subsequent increase in levetiracetam clearance.

#### *3.3. Simulations*

To determine the effects of renal function (expressed by the significant covariate CKD-EPI eGFR) on levetiracetam levels, simulations were performed based on the developed pharmacokinetic model. The utilized model parameters were those listed in Table 2, while two physiological/pathological situations were further simulated: Subjects with normal renal function [23] and subjects with 50% clearance restriction. It can be seen from Figure 4 that increased levetiracetam concentrations are observed with increasing renal function impairment. In healthy subjects, the mean maximum concentration was 31.1 mg/L (Figure 4a), which increased to 37.67 mg/L and 44.12 mg/L in patients with a 50% reduction in eGFR (Figure 4b) and patients with identical physiological characteristics to those in the study (Figure 4c), respectively. This means that the dosing regimens administered resulted in concentrations above the upper value of the reference range of 12–46 mg/L in almost half of the study patients [14].

**Figure 2.** Goodness-of-fit plots for the final best model. (**a**) Observed vs. population predicted by the model concentrations of levetiracetam. The closed circles refer to the (predicted, observed) pairs and the solid line expresses the ideal situation of unity (i.e., y = x). (**b**) Normalized prediction distribution errors (NPDE) vs. time, and (**c**) NPDE vs. concentration. The dotted lines refer to the predicted median (at y = 0) and the 90% predicted percentiles, while the band indicates the 90% prediction interval.

**Figure 3.** Relationship between levetiracetam clearance and the estimated glomerular filtration rate (eGFR) according to the Chronic Kidney Disease Epidemiology collaboration mathematical formula (CKD-EPI eGFR). The eGFR axis is expressed in dual units: mL/min and L/h. The correlation coefficient between levetiracetam clearance and eGFR was equal to 0.76.

**Figure 4.** Simulated concentration vs. time profiles of levetiracetam in normal subjects (**a**), renal impairment 50% (**b**), and for the patients of the study (**c**). In each group, 500 individuals were generated. In the case of the normal renal function subjects (**a**) the utilized model parameters were obtained from study [23]. For the 50% renal impairment the model parameters were similar to (**a**), except for clearance, which was set at 50% of the reported value. For the study patients, the parameter values listed in Table 2 were used. The dosing regimen for each subject was identical to that actually used in the ICU (Table 1).

#### **4. Discussion**

The aim of this study was to investigate the pharmacokinetics of levetiracetam in critically ill intensive care adult patients and to identify possible pathophysiological factors affecting the kinetics of levetiracetam and, consequently, clinical outcome. A clinical study was conducted on 14 patients from an intensive care unit of a tertiary care hospital in whom levetiracetam was administered for either prophylactic or therapeutic reasons. Eleven patients had stable renal function according to the physicians' assessment, while the remaining three had acute renal injury. Having stable renal function does not mean that their renal function was normal, i.e., a glomerular filtration rate greater than 90 mL/min. In these 11 patients, the median eGFR was 97.7 mL/min and three patients (3 of 11) had an eGFR of less than 70 mL/min. In the patient group with acute kidney injury, the median eGFR was 45.3 mL/min. Overall, only five patients had an eGFR greater than 90 mL/min

(i.e., normal renal function), and the remaining nine patients ranged from mild to severe loss of renal function.

Using measurements of serum concentrations in these patients, a population pharmacokinetic model for levetiracetam was developed. The final model for levetiracetam thus derived, referred to a linear one-compartment model with intravenous infusion administration and first-order elimination. This one-compartment model was consistent with other published levetiracetam studies conducted in children, adults, and elderly patients (Table S1) [23–27,41–43]. It should be mentioned that loading doses of levetiracetam were used. However, in the model development and the simulations (Figure 4) performed in this study, loading doses were not considered because they do not alter steady-state plasma levels, but only shorten the time required to reach these plateau levels.

Among the covariates tested, only measures of renal function showed a significant effect on clearance. Among the different eGFR measures explored, CKD-EPI eGFR was selected because it resulted in the lowest *p*-value (*p* = 0.00164, Table 2) and is one of the most commonly used eGFR measures. No significant relationships were found for other scores describing disease severity (such as APACHE II or SOFA), likely due to the small number of patients participating in the study. Levetiracetam clearance in this study (1.45 L/h, Table 2) was lower than the values reported in the literature (from 2.17 L/h to 6.87 L/h, Table S1). This finding may be attributed to the fact that our study included patients with impaired renal function. Moreover, the clearance estimates would also increase with increasing eGFR of our patients, as beta\_eGFR had a positive sign (equal to 0.25, Table 2) [44]. In the study by Sime et al., which also recruited critically ill adult patients, the estimated clearance was 2.51 L/h [7]. Higher levetiracetam clearance estimates (e.g., 5.9 L/h) are reported in studies of epilepsy patients who are not critically ill, such as the large study by Ito et al. [41]. In that study the median eGFR estimate was 97.6 (mL/min)/1.73m2, while in our study the median eGFR estimate of all participants was much lower 81.62 (mL/min)/1.73 m2. The range of eGFR values in our study was from 18.8 mL/min to 123.3 mL/min, whereas for the study of Ito et al. the reported eGFR was as high as 189 mL/min [41]. The issue that our patients had poorer renal function compared with other studies reporting higher levetiracetam clearance studies, may explain the low clearance values found in our study.

The individual fitting plots for the 14 patients are shown in Figure S2. These fitting results demonstrate the adequacy of the descriptiveness of the developed model. Using only three blood samples per patient and data from only 14 subjects, it was possible to develop a model that describes the pharmacokinetics of levetiracetam. For comparison purposes, simulations were performed for the other relevant models in the literature. The model estimates (structural model, mean model parameter, between-subject variabilities, error model) were those reported in the literature (summarized in Table S1) and were kept fixed, while the covariates were related to our study patients. Then, it was attempted to examine the adequacy of fitting of these literature models to the concentration–time data of the 14 patients in this study. VPCs and normalized prediction distribution errors (NPDEs) plots were constructed (Figures S3 and S4 in the Supplementary Material). Figures S3 and S4 reveal that some literature models adequately fit the C-t data of our study, The best overall performance was observed for the models of Karatza et al., Pigeolet et al., Ito et al., and Chhun et al. [23,26,41,43]. However, for all these models, the performance was not better than for our model.

The positive correlation (correlation coefficient equal to 0.76) between eGFR and levetiracetam clearance, shown in Figure 3, indicates a lower elimination capacity in patients with renal dysfunction. This result is consistent with the prescribing information for levetiracetam, which recommends adjusting the levetiracetam dose based on renal function [24,41,42]. At this point, it should be mentioned that several other covariates related to patient disease severity have been tested for their potential role in influencing levetiracetam pharmacokinetics. However, none of these covariates were found to have a significant effect on any of the model parameters of levetiracetam (i.e., *p* > 0.05).

In a further step, the developed population pharmacokinetic model was used to evaluate the effects of impaired renal function on the kinetics of levetiracetam (Figure 4). In the simulated "healthy" subjects, the dosing regimens resulted in levetiracetam levels that were within the therapeutic range (Figure 4a), which is consistent with the literature [40]. The same levetiracetam doses resulted in significantly increased concentration levels in patients with impaired renal function, such as a 50% decrease in clearance and in the study patients (Figure 4b,c), as well as a higher risk of drug toxicity. It is worth noting that this observation is consistent with recommendations in the literature regarding treatment with levetiracetam. It is of note that if seizures cannot be controlled even with high concentrations of levetiracetam, increasing the dose above 2000 mg/day will not benefit the patient; instead, clinicians should consider switching to an antiepileptic drug with a different mechanism of action [23,42]. These findings highlight the importance of dose adjustment in patients depending on their eGFR [45,46]. In general, monitoring levetiracetam serum concentrations as part of therapeutic drug monitoring, ensures more effective and safer dosing regimens.

The relationship between patients' epileptic activity (group 1: no epileptic activity, group 2: possible epileptic activity, group 3: confirmed epileptic activity), APACHE II score, and eGFR levels has also been investigated. Despite the fact that only a few subjects were included in the study, statistically significant (*p* = 0.034) differences in the APACHE II score between the three groups of patients (no epileptic activity, possible, and confirmed epileptic activity) were found. The same groups also differed statistically significantly in terms of Cockcroft–Gault eGFR (*p* = 0.031).

The study had some limitations, one of which was the small number of patients who participated. Due to the small sample size, correlations between variables and their influence on pharmacokinetics may not have been detected (e.g., body weight) [45,46]. Other studies have been published with large numbers of subjects of epileptic adult patients receiving levetiracetam [25,26,40,41]. However, in this study special emphasis was based on exploring the most appropriate index for renal function in the case of critically ill patients. For this reason, this study explored seven indices for glomerular filtration rate (urea, creatinine, renal impairment stage, Cockroft–Gault\_eGFR, MDRD\_eGFR, CKD-EPI\_eGFR, and Jeliffe\_eGFR). Even using the data from this small sample, one covariate (the eGFR values collected from CKD-EPI) was found to perform best compared with all other indices and to play a significant role in the pharmacokinetics of levetiracetam. To increase the convergence of the model, five Markov chains were run in parallel. All available data were replicated, resulting in a number of 50 patients and improving the precision of the estimates.

#### **5. Conclusions**

A clinical trial was conducted on fourteen critically ill patients in an intensive care unit. Levetiracetam was administered intravenously either for prophylaxis or treatment of seizures. Blood samples were collected and plasma levels of levetiracetam were quantified. A population pharmacokinetic analysis was performed that included a number of patientrelated somatometric characteristics and pathophysiological variables. A one-compartment model with intravenous infusion and first-order elimination kinetics was found to best describe the concentration–time data. Several renal function measures were estimated and evaluated in the modeling. Among them, CKD-EPI eGFR showed a statistically significant effect on levetiracetam clearance. With impaired renal function (especially at lower CKD-EPI eGFR values), levetiracetam clearance decreased. Simulations of dosing regimens revealed that even typically administered doses of levetiracetam may result in significantly increased concentrations and a higher risk of drug toxicity in patients with impaired renal function. Therefore, renal function should be considered when developing levetiracetam dosing regimens. Other measures describing disease severity, such as SOFA and APACHE II, had no effect on levetiracetam kinetics. Finally, when analyzing the association between

APACHE II score and epileptic activity, it was found that patients with the lowest APACHE II score had a lower risk of experiencing epileptic seizures.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/app12031208/s1, Table S1: Previously published population pharmacokinetic modeling results for levetiracetam, Figure S1: Goodness-of-fit plots for the final best model (a. Observed vs. predicted by the model individual concentrations of levetiracetam, b. Individual Weighted Residuals (IWRES) vs. time, and c. Individual Weighted Residuals vs. concentration), Figure S2: Levetiracetam concentration vs. time plots for the 14 patients of the study, Figure S3: Visual predictive check plots for the literature models, Figure S4: Normalized prediction distribution errors (NPDE) vs. time and concentration for the literature models.

**Author Contributions:** Conceptualization, S.-L.M., N.M., A.K., D.K., H.S. and V.K.; methodology, S.-L.M., N.M. and V.K.; data curation, R.K. and V.K.; investigation, N.M., A.K., D.K., H.S.; writing original draft, S.-L.M., R.K. and V.K.; supervision, S.-L.M., N.M. and V.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Ethics Committee of Thriassio General Hospital of Elefsina (27 July 2020).

**Informed Consent Statement:** Written informed consent was obtained from the patients to publish this paper.

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

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

#### **References**

### *Article* **Combining Yoga Exercise with Rehabilitation Improves Balance and Depression in Patients with Chronic Stroke: A Controlled Trial**

**Yen-Ting Lai 1,2,3, Chien-Hung Lin 4,5,6, City C. Hsieh 7, Jung-Cheng Yang 2,3, Han-Hsing Tsou 8,9,10, Chih-Ching Lin 4,11, Szu-Yuan Li 4,11, Hsiang-Lin Chan 12,† and Wen-Sheng Liu 4,6,8,9,13,14,\*,†**

	- **\*** Correspondence<EMAIL_ADDRESS>Tel.: +886-2-979305704
	- † These authors contributed equally to this manuscript.

**Abstract:** Background: We combined yoga with standard stroke rehabilitation and compared it to the rehabilitation alone for depression and balance in patients. Methods: Forty patients aged from 30 to 80 who had suffered a stroke 90 or more days previously were divided evenly with age stratification and patients' will (hence not randomized). In the intervention group 16 completed 8-week stroke rehabilitation combined with 1 h of yoga twice weekly. Another 19 patients completed the standard rehabilitation as the control group. Results: The yoga group showed significant improvement in depression (Taiwanese Depression Questionnaire, *p* = 0.002) and balance (Berg Balance Scale, *p* < 0.001). However, the control group showed improvement only in balance (*p* = 0.001) but not in depression (*p* = 0.181). Further analysis showed both sexes benefitted in depression, but men had a greater improvement in balance than women. Depression in left-brain lesion patients improved more significantly than in those with right-brain lesion, whereas balance improved equally despite lesion site. For patients under or above the age of 60, depression and balance both significantly improved after rehabilitation. Older age is significantly related to poor balance but not depression. Conclusions: Combining yoga with rehabilitation has the potential to improve depression and balance. Factors related to sex, brain lesion site and age may influence the differences.

**Keywords:** stroke; yoga; balance disability; depression

**Citation:** Lai, Y.-T.; Lin, C.-H.; Hsieh, C.C.; Yang, J.-C.; Tsou, H.-H.; Lin, C.-C.; Li, S.-Y.; Chan, H.-L.; Liu, W.-S. Combining Yoga Exercise with Rehabilitation Improves Balance and Depression in Patients with Chronic Stroke: A Controlled Trial. *Appl. Sci.* **2022**, *12*, 922. https://doi.org/ 10.3390/app12020922

Academic Editors: Francesco Cappello and Magdalena Gorska-Ponikowska

Received: 1 December 2021 Accepted: 10 January 2022 Published: 17 January 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/).

#### **1. Introduction**

Stroke is highly prevalent in high-income countries (1015–1184 cases per 100,000 people in 2013) [1] Imbalance and post-stroke depression are common among stroke patients and may severely impair their activities of daily living [2,3]. The majority of stroke patients (around 83%) had a balance disability [4]. In a meta-analysis, the percentage of post-stroke depression was about 31% [5], and depression is associated with falls and negatively correlated with functional status [6]. The risk of falls for stroke patients is 1.77 times compared to healthy adults so post-stroke rehabilitation is essential [7].

Exercise increases serotonin which works against depression [8]. Carek et al. reported both aerobic and anaerobic exercises are effective [9]. Increased exercise frequency and group therapy seem be beneficial [10]. The practice of exercises promotes improvements in the levels of depression in people who suffered an ischemic stroke [11]. Franklin et al. reported yoga practices might be also protective [12]. Yoga is widely integrated into physical therapies [13]. However, few reports have focused on the impact of yoga in stroke patients [14–16]. Yoga can lower blood pressure, serum cholesterol and increase in internal reasons for exercise [17,18]. Furthermore, yoga consists of dynamic/static postures to stretch muscles and joints, and improve balance and coordination [19]. Bastille and Gill-Body [14] reported an 8-week intervention (1.5-h yoga sessions two times per week) for patients with chronic post-stroke hemiparesis improved balance. Chronic stroke is defined as the condition 90 days after stroke [20]. About 95% of patients achieved their best neurological recovery within 11 weeks, while minor stroke related to faster recovery [21].

Previous studies enrolled patients who began yoga "after" completing the post stroke rehabilitation program. However, no study focused combining yoga into a standard rehabilitation. We compared depression and balance among chronic stroke patients who underwent combining yoga and standard rehabilitation and those with standard rehabilitation alone.

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

#### *2.1. Ethic Statement*

Institutional Review Board (IRB)/Ethics Committee approval was obtained before the trial began, and the study was conducted in full compliance with the Declaration of Helsinki. The study was approved by the Ethics Committee of National Taiwan University Hospital Hsin-Chu Branch, 201210006IRB. All written consents were obtained from study participants. Clinical Trial Registration Information: Yoga Exercise for Improving Balance in Patients with Subacute and Chronic Stroke, NCT01806922. URL: https://clinicaltrials. gov/ct2/show/NCT01806922 (accessed on 11 February 2014).

#### *2.2. Inclusion and Exclusion Criteria*

The inclusion criteria: (1) Patients with stroke more then 90 days, (2) can stand independently for 1 min, (3) aged from 30 to 80 years. Exclusion criteria: (1) Those under other therapy (such as acupuncture), (2) diagnosis of other psychiatric disease or use of a psychotropic agent, (3) difficult to follow instructions, and (4) other contraindications like cardiopulmonary disease.

#### *2.3. Protocol*

Between onset and enrollment, all patients received the same regular rehabilitation at our hospital. The participants were first stratified into several age groups (5 years apart), then divided each age group by half into the following two groups (first by their will and then by assignment). The intervention group was trained with 60 min of yoga exercise, two times per week, combined with standard stroke rehabilitation for 8 weeks. The control group only underwent standard rehabilitation (Figure 1).

**Figure 1.** Flow diagram of the study.

Yoga Exercise Combined Standardized Rehabilitation Program (Intervention Group): In the intervention group, yoga was taught by a registered yoga therapist. A standard protocol was developed with 10 min of warm-up, 40 min of main activity, and 10 min cool-down exercises. All sessions focused on breathing, deep relaxation, mediation, posture, and trunk/limb stretching based on a study showed yoga improved balance in chronic stroke patients [15] (see Supplementary Materials Table S1–S3).

Standardized Rehabilitation Program (Control Group):All participants participated in the program based on the Taiwan Guideline for Stroke Rehabilitation [22] focused on dynamic balance exercises and weight-shift training given by the same therapists.

#### *2.4. Assessment of Outcome*

The outcome (depression and balance) were measured with the Taiwanese Depression Questionnaire (TDQ) [23] and Berg Balance Scale (BBS) [24] by the same therapist. (For TDQ: the lower the score, the better the mood. As for BBS: the higher the score, the better the balance.)

The TDQ is an 18-item questionnaire, and each item applies a 4-point scale for the response. It is a culturally specific self-rating instrument for screening of depression. Subjects are guided to rate each item on a scale from 0 to 3 based on "how often you felt during the past week". TDQ scores range from 0 to 54. Previous studies indicate good psychometric properties in the TDQ [25,26]. The cutoff scores for the screening of depression was 15 [27].

Balance was assessed with the Berg Balance Scale (BBS), which is the gold standard for functional balance testing [28]. It is a 14-item scale and can measure static and dynamic balance after stroke with good reliability and validity [28]. The total score is 56 points. The risk of falls in older adults was high when the score was less than 36 [29], and the

risk increased by 6% to 8% with every one-point reduction within the range of 46 to 54 points [30].

#### *2.5. Statistical Analysis*

We performed statistical analyses using SPSS for Windows (Version 19.0, IBM Corp, Armonk, NY, USA). Continuous data are expressed as mean ± standard deviation. The independent *t*-test was used to compare the differences between the two groups (indicated with dotted arrows in the figure). We analyzed the measured values before and after the treatment by paired *t*-tests (indicated with bold arrows in the figure). A *p*-value of less than 0.05 is considered to be statistical significant (indicated with an asterisk in the figure). The sample size was determined based on an effective size to detect differences in different groups. If we permitted a 5% chance of a type I error (α = 0.05), with a power of 80%, assuming the difference among the two groups was at least equal to the standard derivation, then approximately 32 patients would be required, while 35 patients completed the study. Subgroup analysis was performed on sex, brain lesion site and age.

#### **3. Results**

#### *3.1. Primary Analysis of the Effect for Adding Yoga*

Forty people enrolled initially with 20 to the intervention (yoga) group and 20 to the control group (Figure 1). In the yoga group, 16 completed the study. (One failed due to feeling exhausted, and three due to loss of follow-up thus being excluded). In the control group, 19 completed the study, and one loss of follow-up due to recurrent stroke. The completion rate of did not differ significantly. Age, sex, and brain lesion site (Table 1) were also similar. The TDQ score did not differ at baseline, whereas BBS scores were significantly better in the yoga group (Table 1). Overall, the TDQ of all patients decreased (improved) (15.23 ± 9.91 to 11.71 ± 8.66; *p* = 0.001) and the BBS increased (improved) from 44.71 ± 7.64 to 46.60 ± 7.46 (*p* < 0.001) after intervention.

**Table 1.** Comparison of baseline clinical and 8-week post-assessment of depression (Taiwanese Depression Questionnaire, TDQ) and balance (Berg Balance Scale, BBS), between control and experimental groups.

Data are given as means ± SD unless otherwise indicated. † One patient was eliminated in this analysis due to this patient suffered from bilateral stroke. \* *p* < 0.05, T onset: days after stroke. D0: for baseline TDQ; D1 for post-assessment TDQ; delta D: the change from D0 to D1. B0: for baseline BBS; B1 for post-assessment BBS; delta B: the change from D0 to D1.

The TDQ scores improved significantly only in the yoga group (14.06 ± 8.85 to 8.25 ± 5.19, *p* = 0.002). In contrast, patients in the control group had little improvement (*p* = 0.181) (Figure 2A) (bold black arrow). The baseline TDQ did not differ (*p* = 0.531). But the post-study TDQ was significantly lower (8.25 ± 5.19 vs. 14.63 ± 9.99, *p* = 0.028) in the yoga group (Figure 2A) (Dotted arrow).

**Figure 2.** (**A**): Bold black arrows show TDQ score significant dropped (improved) in the yoga group after the intervention compared with their own baseline (*p* = 0.002). Dotted arrow showed the post intervention TDQ score is significantly lower in the yoga group compared to the control standard rehabilitation group. (*p* = 0.028). (**B**): Bold black arrows show BBS score significant increased (improved) in yoga (*p* < 0.001) and control group (*p* = 0.001) after the intervention compared with their own baseline BBS score. The dotted arrow shows the pre- and post-study values of the BBS in the yoga group are higher than that of the control group (both *p* < 0.001) Screwed dotted arrow shows the post intervention increase of BBS score is significantly larger in the yoga group compared to the control standard rehabilitation group (*p* = 0.029). (\*: *p* < 0.05).

Both groups showed improvement in balance. The BBS score increased from 48.56 ± 4.42 to 51.13 ± 4.25 (*p* < 0.001) in the yoga group, while from 41.47 ± 8.01 to 42.79 ± 7.85 (*p* = 0.001) (Figure 2B) for the control group. (Bold black arrows) Although the pre- and post-study BBS scores in the yoga group were higher (dotted arrows), the improvement in the yoga group is also significantly larger than that of the control group (2.56 ± 1.96 vs. 1.31 ± 1.45, *p* = 0.029) (Table 1) (Figure 2B) (screwed dotted arrow).

#### *3.2. Subgroup Analysis*

The experiment group and control group were combined and divided by gender (female vs. male), brain lesion side (left vs. right) and age (age <60 vs. 60 or more) sequentially to see how these factors affected the result of rehabilitation (Table 2).

**Table 2.** Subgroup analysis for sex, lesion side, and age before and after the rehabilitation in comparison Taiwanese Depression Questionnaire (TDQ) scores and Berg Balance Scale (BBS) scores.

\* *p* < 0.05, T onset: days after stroke. D0: for baseline TDQ; D1 for post-assessment TDQ; delta D: the change from D0 to D1. B0: for baseline BBS; B1 for post-assessment BBS; delta B: the change from D0 to D1.

#### 3.2.1. Gender

The mean age and the duration since stroke did not statistically differ by sex (Table 1). The TDQ scores improved after the study in both men and women. The TDQ scores dropped (improved) from 16.73 ± 9.18 to 12.93 ± 8.10 (*p* = 0.020) for women and from 14.10 ± 10.51 to 8.10 ± 9.16 (*p* = 0.027) for men (Figure 3A) (bold black arrow). For balance, the BBS score among all women improved from 44.60 ± 8.16 to 45.73 ± 8.53 (*p* = 0.009) and from 44.80 ± 7.08 to 47.25 ± 7.05 among men (*p* < 0.001) (Figure 3B). Thus, rehabilitation was effective despite sex. However, significantly greater improvement in BBS was observed in men (*p* = 0.029) (Figure 3B) (screwed dotted arrow). Balance among men seems to benefit more than that of women.

#### 3.2.2. Brain Lesion Side

Further analysis of the stroke site in the brain is shown in Table 2 and illustrated in Figure 3B. For depression, patients with left-sided brain lesions showed significant improvement (TDQ scores decreased from 13.47 ± 9.28 to 8.94 ± 5.68, *p* = 0.006) (Figure 3C) (bold black arrow). In contrast, patients with right-sided brain lesions showed no significant improvement in TDQ (*p* = 0.076) (Figure 3C).

For balance, unlike depression, patients showed significant improvement despite the lesion site. Overall, BBS scores of patients with left-brain and right-brain lesions increased from 46.41 ± 6.21 to 48.82 ± 6.24 (*p* < 0.001) and from 44.12 ± 7.34 to 45.35 ± 7.80 (*p* = 0.003), respectively (Figure 3D) (bold black arrow). In summary, TDQ scores for depression improved only in patients with left brain lesions, but balance improved in all subgroups despite the lesion site.

**Figure 3.** (**A**): Bold black arrow shows TDQ score significant dropped (improved) in the male and female groups after the intervention compared with their baseline TDQ score. (**B**): The improvement of TDQ is mainly contributed by the yoga groups in both genders. (**C**): Bold black arrow shows TDQ score significant dropped (improved) in the left brain lesion group after the intervention compared with their baseline. (*p* = 0.006) (**D**): The improvement of TDQ is mainly contributed by patients with left brain lesion in the yoga groups. (*p* = 0.011) (**E**): Bold black arrow shows TDQ score significant dropped (improved) in the young group (*p* = 0.032) and old group (*p* = 0.015) after the intervention compared with their baseline. (**F**): The improvement of TDQ is mainly contributed by young patients in the yoga groups (*p* = 0.010). (\*: *p* < 0.05).

#### 3.2.3. Age

Table 2 also shows the age effect. We divided the patients by the age of 60 because the mean age is about 60 years (55.79 ± 7.81). Overall, depression improved in both the young and older groups. The TDQ score dropped from 15.74 ± 8.93 to 12.32 ± 7.97 (*p* = 0.032) in the young and from 14.63 ± 11.23 to 11.00 ± 9.64 (*p* = 0.015) in the older group (Figure 3E) (bold black arrow).

As for balance, both young and older patients benefitted from rehabilitation. Generally, the BBS score improved from 48.58 ± 4.77 to 50.3 ± 4.65 (*p* = 0.001) in young patients and from 40.13 ± 7.55 to 42.13 ± 8.18 (*p* < 0.001) among older patients (Bold black arrows). Additionally, the BBS scores of older patients were much lower than those of the younger ones. (Figure 3F) (dotted arrows).

#### *3.3. Factors Correlated with Balance*

Further analysis showed age is a big contributor to poor balance (linear regression for age vs. baseline BBS score, *p* = 0.005, R2 = 0.464; age vs. post-intervention BBS score, *p* = 0.006, R2 = 0.456) (Figure 4).

**Figure 4.** (**A**) Age is associated with lower baseline balance score (BBS) (*p* = 0.005, R2 = 0.464), (**B**) Age is associated with lower post intervention balance score (BBS) (*p* = 0.006, R2 = 0.456) (\*: *p* < 0.05)).

Although the baseline TDQ and BBS scores are not correlated (*p* = 0.25), the postintervention TDQ score is significantly correlated to the BBS score (*p* = 0.016) (Figure 5).

**Figure 5.** (**A**) Baseline balance score is not associated with baseline depression score (TDQ) (*p* = 0.250, R2 = 0.040), (**B**) post intervention balance is associated with post intervention depression score (TDQ) (*p* = 0.016, R2 = 0.163). (\*: *p* < 0.05).

#### **4. Discussion**

All patients had significant improvement in depression and balance after rehabilitation (Table 1). Overall, a decrease in TDQ score (*p* = 0.001) and an increase in BBS score (*p* < 0.001) were noted. We found the TDQ scores of patients in the yoga group were much lower than those in the control group after the intervention while the baseline TDQ scores were similar (Figure 2A). Only patients in the yoga group showed significant improvement in depression after the intervention which indicates yoga provided anti-depression effect.

Because this was not a randomized trial, the baseline BBS score is much higher in the yoga group (with a better physical condition) compared to the control group even before intervention. BBS scores remained much higher after the intervention (Figure 2B). The reason may be that healthier people have higher intention to be volunteers for epidemiological research or a challenging task and have better survival rates than the general population, thereby producing bias referred to as the "healthy volunteer effect" [31]. Although both groups showed improvements in BBS scores (Figure 2B), the change of improvement in BBS was also significantly larger in the yoga group (*p* = 0.029). This finding indicates yoga still provided extra benefits compared to standard rehabilitation.

We compared the differences in the effects of rehabilitation on both sexes (Table 2 and Figure 3). Overall, we found that depression was significantly improved equally (Figure 3A). Unlike depression, balance showed a significantly greater improvement in men generally than in women (*p* = 0.029) (Figure 3B). One possible explanation is that men have poorer flexibility than women due to a lack of hormones such as estrogen or relaxin. Relaxin is known for relaxing the tendons and joints in pregnancy and delivery in women [32]. Yoga has been proven to improve flexibility [33], and flexibility is related to balance [34]. It may be that women have better flexibility in the first place, so the improvement of flexibility and balance after yoga may not be as good as man. In other words, the mind factor (anti-depression) of yoga is similar for both sexes but the body factor (balance) is quite different with better improvements in men.

The effect of brain lesion site is shown in Table 2 and Figure 3B. Only patients with left-sided brain lesions showed a significant improvement in depression (Figure 3C). This is probably due to left-sided brain lesion being associated with aphasia because the linguistic center is located the left brain while aphasia is highly associated with depression [35]. Yoga may lead to an improvement in aphasia [36]. Unlike depression, balance was uniformly improved in all subgroups of patients, despite of brain lesion side. These results are compatible with the brain structure whereby motor cortex is symmetrical, but the language center is located in the left brain [35].

Age is also an influencing factor. For patients younger than 60, depression (*p* = 0.032) and balance (*p* = 0.001) significantly improved. (Figure 3E). For older patients aged 60 or more, depression showed significant improvement (*p* = 0.015), and balance improved as well (*p* < 0.001) (Figure 3F). Therefore, rehabilitation is important for both young and older patients.

Further analysis of factors associated with balance, we found the BBS score is much lower in the older group, as age is significantly inverse related to the balance (Figure 4). This could be because young patients have better neuroplasticity and muscle power [37]. Aging is a big contributor to poor balance (linear regression for age vs. baseline BBS score, *p* = 0.005, R2 = 0.464; age vs. post-intervention BBS score, *p* = 0.006, R2 = 0.456) (Figure 4A,B). Nonetheless, age is not associated with depression (Age vs. baseline TDQ score, *p* = 0.551; age vs. post-intervention BBS score, *p* = 0.439). Therefore, old age may not guarantee depression.

On the other hand, although the baseline depression score (TDQ) and balance (BBS) scores are not correlated (*p* = 0.25), the post-intervention TDQ score is significantly correlated to the BBS score (*p* = 0.016). This may indicates that better post-intervention physical condition may contribute to a better quality of life and mood after completing the rehabilitation (Figure 5A,B). However, there are more factors other than age (such as gender, brain lesion side) affecting the outcome of rehabilitation.

In summary, for the subgroup analysis of overall rehabilitation effect, we found both sexes benefited in depression and men benefited more in balance than women. Patients with left-brain lesions improved more in depression than those with right lesions after rehabilitation, while balance equally improved regardless of the lesion site. For patients under or above the age of 60, depression and balance both significantly improved after rehabilitation. Older age was significantly related to poor balance but not depression. Post-intervention higher balance (BBS) score was significantly correlated to lower TDQ score and thus a better mood.

Our first limitation was that this was a pilot study with small recruitment number, but we still obtained many significant results. Second, the inclusion criteria limit each subject must be able to stand independently for 1 min, which indicated participants were all relatively mildly functionally impaired. Whether the results can be applied to more disabled people is unknown. Third, in order to respect participants' will and their capability, our study was not randomized. (However, if done in a randomized manner, the completion rate could be much lower in the intervention group.) Also, the long-term effect after 8 weeks was not confirmed. A larger randomized controlled study with longer follow up period is necessary in the future. Finally, we should add one more group with standard rehabilitation of 300 min (5 h) each week to see which group performs better. However, patients always get bored and become impatient during repetitive exercise in the real clinical setting. That is why we add recreational sport (such as yoga) as an adjuvant exercise because yoga is relaxing and interesting.

Because of the COVID-19 pandemic, group exercise may increase the infection risk. Whether remote video-assisted yoga training provides the same benefit as the current group exercise requires further research.

#### **5. Conclusions**

Generally, improvements were noted for balance in patients with standard rehabilitation but improvement in depression was only observed among patients when yoga was combined. Yoga provided more benefit for men than women in balance, for patients with left side brain lesion in depression, and for patients aged under or above 60 in both balance and depression. Older age was significantly related to poor balance but not depression. Post-intervention higher balance (BBS) score was significantly correlated to lower TDQ score and thus a better mood. Our study showed that combining yoga in the current standard rehabilitation is safe and has the potential to improve mood and balance, which provides useful information for stroke patients in tailor-made rehabilitation combing yoga.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/app12020922/s1. Table S1: Standard stroke rehabilitation; Table S2: Yoga class schedule (16 yoga sessions) in each week; Table S3: Yoga pose.

**Author Contributions:** Y.-T.L. made substantial contributions to the conception and design; acquisition, analysis and interpretation of data; drafting of the manuscript. C.-H.L., H.-H.T., S.-Y.L. and C.-C.L. revised the manuscript critically for important intellectual content. C.C.H. and J.-C.Y. provided clinical care for the patient and critically revised the manuscript for important intellectual content. W.-S.L. and H.-L.C. served as primary consultants in the management of the patient and made substantial contributions to the conception and design; analysis and interpretation of data; drafting of the manuscript; and revising the manuscript critically for important intellectual content; and gave final approval of the version to be published. All authors have read and agreed to the published version of the manuscript.

**Funding:** This trial was supported by the Department of Physical Medicine and Rehabilitation, National Taiwan University Hospital Hsin-Chu Branch (Protocol No.: HCH-102-22), the Ministry of Science and Technology (MOST 103-2410-H-134-022-MY), Taiwan and Taipei City Hospital (TPCH-108-27), and Department of Health of Taipei City Government (10901-62-023).

**Institutional Review Board Statement:** Institutional Review Board (IRB)/Ethics Committee approval was obtained before the trial began, and the study was conducted in full compliance with

the Declaration of Helsinki. The study was approved by the Ethics Committee of National Taiwan University Hospital Hsin-Chu Branch, 201210006IRB. Clinical Trial Registration Information: Yoga Exercise for Improving Balance in Patients with Subacute & Chronic Stroke, NCT01806922. URL: https://clinicaltrials.gov/ct2/show/NCT01806922 (accessed on 11 February 2014).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study. Written informed consent has been obtained from the patient(s) to publish this paper.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

### *Article* **Development of a Simple Spheroid Production Method Using Fluoropolymers with Reduced Chemical and Physical Damage**

**Hidetaka Togo 1, Kento Yoshikawa-Terada 2, Yudai Hirose 1, Hideo Nakagawa 1, Hiroki Takeuchi 2,\* and Masanobu Kusunoki 1,\***

**Abstract:** Establishing an in vitro–based cell culture system that can realistically simulate in vivo cell dynamics is desirable. It is thus necessary to develop a method for producing a large amount of cell aggregates (i.e., spheroids) that are uniform in size and quality. Various methods have been proposed for the preparation of spheroids; however, none of them satisfy all requirements, such as cost, size uniformity, and throughput. Herein, we successfully developed a new cell culture method by combining fluoropolymers and dot patterned extracellular matrix substrates to achieve size-controlled spheroids. First, the spheroids were spontaneously formed by culturing them twodimensionally, after which the cells were detached with a weak liquid flow and cultured in suspension without enzyme treatment. Stable quality spheroids were easily produced, and it is expected that the introduction and running costs of the technique will be low; therefore, this method shows potential for application in the field of regenerative medicine.

**Keywords:** spheroid; organoid; van der Waals force; fluoropolymers; adhesion; microstamp

#### **1. Introduction**

Cell culture, an important technique in basic and preclinical research, is a cost-effective way of reducing the number of experimental animals during the drug discovery process [1]. Although two-dimensional (2D) cell culture is a universal and valuable method, it has its limitations [2]; for example, the cells proliferate and grow in a single layer and thus lack the cell–cell and extracellular matrix (ECM) interactions present in native tissues and cells. In addition, the cells are stretched and undergo a cytoskeletal rearrangement, acquiring artificial polarity and expressing abnormal genes and proteins [3]. Unlike tissues in an in vivo environment, 2D systems cannot provide a complex and dynamic microenvironment for cells, leading to potentially invalid findings [4,5]. For these reasons, the majority of compounds discovered in 2D cell culture during drug discovery fail in clinical trials [6]. Therefore, it is important to establish an in vitro cell culture system that can more realistically simulate the dynamics of cells in vivo.

Three-dimensional (3D) culture systems are an excellent in vitro model that can mimic in vivo processes, such as embryogenesis, morphogenesis, and organogenesis [7]. Spheroids and organoids are cellular aggregates with complex cell–cell adhesion that generate gradients of nutrients, gases, growth factors, and signaling factors. These structures can replicate the cellular microenvironment observed in real tissues and simulate in vivo cell dynamics more realistically. Human induced pluripotent stem cells (iPSCs), the discovery of which has helped stem cell research make great strides, can be used to generate spheroids and organoids, which were employed in the elucidation of pathological conditions [8–12]

**Citation:** Togo, H.;

Yoshikawa-Terada, K.; Hirose, Y.; Nakagawa, H.; Takeuchi, H.; Kusunoki, M. Development of a Simple Spheroid Production Method Using Fluoropolymers with Reduced Chemical and Physical Damage. *Appl. Sci.* **2021**, *11*, 10495. https://doi.org/ 10.3390/app112110495

Academic Editor: Francisco Arrebola

Received: 19 September 2021 Accepted: 4 November 2021 Published: 8 November 2021

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

**Copyright:** © 2021 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/).

and drug discovery [8,13–15]. Recently, organoids derived from iPSCs generated from patients with intractable diseases were used to elucidate the disease mechanisms [16,17]. In the field of drug discovery, the validation of drug efficacy using organoids derived from human iPSCs in the early stage of in vitro testing is expected to reduce the possibility of candidate compound failure before initiating clinical trials. Indeed, drug testing using iPSC-derived organoids has shown comparable results to those of drugs already used for treatment [18]. Because the use of spheroids and organoids is expected to increase in the future, developing methods for their efficient production is important.

Various 3D culture methods have been tried and implemented. Conventional 3D culture methods can be broadly classified into four types: micropatterned forced flotation, hanging drop, swirling culture or spinner flask, and force-driven methods [19]. The micropatterned forced flotation method involves well plates in which cells deposited inside the microwells adhere to each other in three dimensions [20–22], as well as low-adhesion plates in which areas with different cell adhesion properties are patterned to prevent cell aggregation [23]. Microfluidic systems are a type of micropatterned forced flotation that combines both the target cells and matrix components to create spheroids on microfluidic chips with microgel beads [24,25]. The hanging drop method utilizes the surface tension of the culture medium to form droplets, while the interface with the air is used as a microwell [26–28]. The spinner flask or swirling culture (bioreactor) is a method of spheroidization via collision and contact between cells through continuous swirling and agitation of the culture medium in which cells are dispersed [29,30]; this method also includes microfluidic types on a smaller scale [31–33]. Force-driven methods include magnetic methods where nanoparticles, such as iron oxide, gold, and poly-L-lysine, are attached to the cell surface. Spheroids are formed via the following: magnetic forces [34]; electric methods in which an external force, such as an electric field, is applied [35]; and acoustic tweezers [36]. Other methods have also been developed in recent years; for example, 3D printing can artificially create 3D shapes [37–39], using a computer-controlled arrangement of heterogeneous cells with a gel-like culture medium as the filler. Although these methods exhibit excellent features, such as dimensional control of the spheroid, reduction in damage and irritation to the cells, simplicity of the procedure, high success rate, mass production, high throughput, and low cost, none of them satisfy all these requirements simultaneously. A new type of classification is microfluidic devices. Microfluidic devices are devices that have microscopic grooves or tubes within. These are used for various applications, such as forming small bioreactors [32,40] or circulatory systems [31] to recapitulate in vivo, using flow paths. Some of these microfluidic methods are included in the micropatterned forced flotation method.

Fluoropolymers have a low surface free energy and are generally regarded as materials that do not allow cell adhesion. They have a high biocompatibility and inertness, due to these properties. Therefore, fluoropolymers are used in artificial blood vessels and other components that are implanted in the body. In general environments, these polymers are very stable, due to the high binding energy of C–F bonds, and they are, thus, a safe and non-toxic material for use in the living body. However, their high inertness prevents cell adhesion, which is a problem for biomaterials. As a result, various surface modifications have been investigated to take advantage of their high biocompatibility and inert properties [41–46]. To improve adhesion, plasma treatment [47–51], UV irradiation [52–55], γ-radiation [56,57], ion introduction [58–61], and polydopamine treatment [62–66] have been employed for surface modification; however, most of these methods require high energy to break the strong C–F bonds and add functional groups to improve adhesion.

We are working on the control of cell adhesion by manipulating the surface properties of fluoropolymers for 3D cell patterning and have found that extracellular matrix (ECM) substrates deposited on fluoropolymers show weak cell adhesion in this study. In this study, we took advantage of this property and developed a new spheroid fabrication method that simultaneously satisfies all of the above requirements for spheroid production by forming an ECM dot pattern on the fluoropolymer, which allows cells to adhere on

it, retain it during cell growth, and detach from the fluoropolymer without employing damaging enzymatic treatments.

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

#### *2.1.  CYTOP™ is a fluoropolymer product with properties similar to polytetrafluoroethylene (PTFE) [67,68]. This resin has fluoropolymer properties, such as water repellency, oil repellency, and chemical resistance, as well as high transparency, which are difficult to obtain with crystalline resins, such as PTFE. It can also be used for cell culture observation using an inverted microscope because it has no fluorescence properties. The material was immersed in a dip solution consisting of CYTOP™ and a 16-fold dilution of a thinner solution (CT-Solv.100; AGC Chemicals) to coat the glass substrate (Table 1). The pull-out and curing conditions were as recommended by the manufacturer. The material was pulled out at a speed of approximately 1 mm/s and cured in an oven at 100 ◦C for at least 90 min.

#### *2.2. Preparation of a Stamp for Patterning Formation*

The dot patterns of the ECM were fabricated using a microstamp made out of polydimethylsiloxane (PDMS; SYLGARD™184; Dow Toray, Tokyo, Japan). PDMS has high biological safety [69–73] and is also used in engineering to transfer fine patterns [74]. To fabricate the microstamp, PDMS mixed with a hardener was poured into a mold made of PTFE, after which the mold was defoamed. After defoaming, the stamps were heat-cured at 90 ◦C for 90 min. The dimensions of the stamps were 100 μm in spheroid diameter (the size at which spheroid necrosis does not occur [75–77]), 800 μm in dot size diameter, and a 7 × 7 matrix arrangement with a pitch of 1500 μm. The dot diameter was defined as the area where the volume of a cell sheet cultured in a monolayer on the dot pattern is equivalent to a spheroid 100 μm in diameter, assuming that one cell is a sphere with a diameter of 10 μm. The dimensions of the molds and microstamps were measured using an OLS-5000 (Olympus, Tokyo, Japan).

#### *2.3. Patterning of the Extracellular Matrix*

ECM substrates, including Matrigel Growth Factor Reduced (354230; Corning, NY, USA), iMatrix-511 silk (892021; Nippi, Tokyo, Japan), fibronectin (33016015; Thermo Fisher Scientific, Waltham, MA, USA), vitronectin (A31804; Thermo Fisher Scientific, Waltham, MA, USA), and collagen IV (ASC-4-104-01; Nippi), were diluted in RPMI 1640 (30264-56; Nacalai Tesque, Kyoto, Japan) 2-, 4-, 8-, and 16-fold. PDMS stamps dipped in ECM diluent were stamped onto the substrates, after which the ECM was transferred to the fluoropolymer surface with the formed pattern and incubated at 37 ◦C and 5% CO2 for 1 h. After incubation, the substrates were washed with RPMI 1640.

#### *2.4. Cell Culture*

HepG2 (TKG02058) and MCF7 (TKG0479) cells were provided by the Cell Resource Center for Biomedical Research, Institute of Development, Aging and Cancer, Tohoku University, Japan. The cells were cultured in Dulbecco's modified Eagle's medium (DMEM; 5919, Nissi, Tokyo, Japan.) supplemented with 10% fetal bovine serum (FBS; BioWest, Nuaillé, France), 2 mM L-Alanyl-L-glutamine solution (01102-82, Nacalai Tesque) and Penicillin-Streptomycin Mixed Solution (26252-94, Nacalai Tesque) (10%FBS/DMEM.) for 37 ◦C at 5% CO2 and 95% air. The cells were passaged every 3 days at a ratio of 1:3, using TE dissociation buffer containing 0.25% trypsin (35547-64; Nacalai Tesque) and 0.04% EDTA (15105-35; Nacalai Tesque).

#### *2.5. Spheroid Preparation*

The ECM was transferred to the fluoropolymer-coated substrate, using a microstamp (Figure 1a). Cells were then evenly seeded at 1.0 × 104 cells/cm2 on the substrate and allowed to evenly adhere (Figure 1b). After 7 days, the cells proliferated and formed colonies only on the ECM-coated area (Figure 1c). These cell colonies were detached by gently pipetting, followed by culturing in suspension for 7 days until spheroid formation (Figure 1d). All cultures were performed on fluoropolymer-coated glass substrate in a 5-well petri dish (NM-CD-5F, Naka Medical, Tachikawa, Japan).

**Figure 1.** Schematic diagram of the spheroid production process. (**a**) Matrix transferred by stamping. (**b**) Cell adhesion occurring one day after seeding. (**c**) Cell colony formation 7 days after seeding. (**d**) After detaching the cells, spheroids formed in suspension by day 14.

Spheroids were produced by the hanging drop method. The average number of cells in the dot pattern colonies was estimated, and they were seeded to form a spheroid. Briefly, the cells were grown until confluent. Single cells were then dissociated using TE buffer. The number of cells was counted, and spheroids were prepared by seeding 8 × 103 cells/drop using the hanging drop method.

#### *2.6. Evaluation of Cell Patterning*

All ECM and cell patterning were observed using a phase-contrast microscope and then photographed and dimensioned using FLOYD-4K (Wraymer, Osaka, Japan). The diameter of the ECM was measured, but as cell colonies and spheroids do not form a perfect circle, the equivalent diameter of a circle was calculated from the area. Photographs were acquired of the coated ECM dimensions, colony size after cell seeding (day 1), cell sheet dimensions left as a pattern on the ECM (day 7), and spheroid size (days 10–14).

#### *2.7. Measurement of Spheroid Viability Using Trypan Blue Staining*

The cells of the 2D culture and spheroids were dissociated using TE buffer, and a single cell suspension was obtained. The cells in the suspension were then stained with trypan blue (35535-02, Nacalai Tesque) to determine the number of viable cells. Student *t*-test was performed to compare survival rates between the 2D culture and spheroids.

#### **3. Results**

#### *3.1. Water Repellency of Fluoropolymers*

The contact angle measurements were conducted to measure the surface free energy, which is mainly responsible for fluoropolymer inertness in cells. The contact angles for a 2-, 4-, 8-, and 16-fold dilution were 113.36◦, 115.00◦, 114.31◦, and 112.58◦, respectively, showing good water repellency, even at a 16-fold dilution. In consideration of economic efficiency, we used a 16-fold dilution of CYTOP™ for coating (Table 1).

**Table 1.** Contact angle measurement results for the different CYTOP™ dilutions.

#### *3.2. Dimensional Reproducibility of the Microstamps*

After fabricating the microstamp with PDMS, we measured the dimensions of the top surface of the dots and found an average of 825.24 μm, σ of 4.96 μm, and the shape diameter that was close to the target of 800 μm with small variations (Figure 2).

**Figure 2.** Observation of the top surface of a micro stamp using a laser microscope. (**a**) Optical imaging of a microstamp; (**b**) 3D image by a laser.

#### *3.3. Culturing HepG2 Cells on the Matrigel Pattern*

Matrigel was transferred to the fluoropolymer-coated substrate using a microstamp with an average pattern diameter of 898.74 ± 140.69 μm (Figure 3a,e: matrix transferred using a microstamp). HepG2 cells were seeded on this substrate and found to be uniformly adhered to the entire surface by day 1 (Figure 3b). The medium was gently changed every 2–3 days to avoid detaching the cells. By day 7, the cells had formed circular colonies based on the stamp pattern with an average colony diameter of 930.11 ± 178.47 μm (Figure 3c,e). After detaching the colonies from the substrate by weak pipetting and culturing in suspension, they spontaneously formed spheroids by day 8 and completely spheroidized by day 14. The average spheroid size was 466.23 ± 75.25 μm (Figure 3d,e, and Supplementary Figure S1).

**Figure 3.** Bright-field image of HepG2 cultured on Matrigel dots and chart of size on every stage. (**a**) Stamped ECM (Matrigel) on the fluoropolymer surface. (**b**) HepG2 cells 1 day after seeding. (**c**) HepG2 cells grown to a colony of the same size as that of the Matrigel pattern on day 7. (**d**) Day 14 spheroids derived from detached HepG2 cells grown in suspension after culturing in 2D for 7 days on Matrigel. (**e**) The average diameter of the dot patterns during each step as shown in (**a**–**d**). Error bars indicate the standard deviation. (**f**) Viability confirmation by trypan blue staining. Student *t*-test was performed to compare survival rates between the 2D and spheroid (*p* = 0.05). *p* value < 0.05 was considered statistically significant. Data are expressed as means ± standard deviations. The asterisks (\*) indicate a statistically significant difference. (**g**) HepG2 on days 1 and 7 following the hanging-drop method.

> Trypan blue staining was performed to determine the viability of the prepared spheroids. The survival rate of the HepG2 cells as estimated by trypan blue staining was 93.95 ± 6.56% in 2D culture and 82.84 ± 6.99% in the spheroids derived by our method (Figure 3f). Next, spheroid formation was monitored for 7 days using the hanging-drop method; however, no HepG2-derived spheroids were observed (Figure 3g).

#### *3.4. Culturing HepG2 Cells on Other Extracellular Matrix Patterns*

We used other ECM substrates and the MCF7 cell line to determine whether spheroids formed using other combinations besides Matrigel and HepG2 cells. Using the other ECM substrates, including iMatrix (fragmented laminin), vitronectin, fibronectin, and collagen IV dot patterns also formed. The seeded HepG2 cells successfully adhered to these ECM substrates, after which colonies formed by day 7 and spheroids by day 14 (Figure 4(a1– a5)). Similar to HeG2 cells, seeded MCF7 cells formed colonies 7 days after seeding and

spheroids by day 14 after the suspension culture (Figure 4(b1–b5)). When the cells were continuously cultured in 2D without detaching by day 7, the cells continued to proliferate and became over confluent on the ECM dots by day 10. Subsequently, the cells no longer proliferated in the monolayer, due to the influence of non-adhesive areas outside the colony pattern but instead proliferated toward the upper part of the colony and partially formed a spheroid. This phenomenon was also observed with HepG2 cells, indicating that it is difficult to control spheroid diameter within the area of the pattern when using cell lines with high proliferation ability.

**Figure 4.** Spheroids derived from HepG2 or MCF7 cells after culturing in 2D for 7 days on different ECM substrates. HepG2 cells seeded on Matrigel (**a1**), Matrix (**a2**), fibronectin (**a3**), vitronectin (**a4**), and collagen IV (**a5**). MCF7 cells seeded on Matrigel (**b1**), Matrix (**b2**), fibronectin (**b3**), vitronectin (**b4**), and collagen IV (**b5**). (**c**) Viability confirmation by trypan blue staining. Student *t*-test was performed to compare the survival rates between the 2D and spheroid (*p* = 0.01). *p* value < 0.05 was considered statistically significant. Data are expressed as means ± standard deviations. The asterisks (\*) indicate a statistically significant difference. (**d**) MCF7 on days 1 and 7 following the hanging-drop method.

> In MCF7 and matrigel, trypan blue staining was performed to determine the viability of spheroids from colonies; the survival rate of MCF7 in trypan blue was 96.09 ± 2.82% in 2D culture and 87.89 ± 3.62% in the spheroids derived by our method (Figure 4c). Next, spheroid formation was monitored for 7 days, using the hanging drop method. MCF7 derived spheroids with a diameter of approximately 100–200 um were observed on day 1; however, no larger-sized spheroids were observed (Figure 4d).

#### **4. Discussion**

In the present study, all tested ECM substrates were successfully scaffolded on a fluoropolymer, and both HepG2 and MCF7 cell lines formed spheroids after culturing, without negatively affecting the cell population. The key to achieving this was the combination of the fluoropolymer and ECM dot pattern. Fluoropolymers are generally known to be difficult regarding adhesion; thus, many studies have attempted to improve adhesion by surface modification [78]. However, we showed that cells can adhere to the substrate via an ECM that can be detached easily by pipetting.

The survival rate of the spheroids derived using our method was significantly lower than that of the cells in the 2D culture. The reason for this was speculated to be the

lack of space for cell growth inside the spheroid and the lack of nutrient circulation in the culture medium. For cancer spheroids, however, a slight decrease in survival are permitted; there are, accordingly, no disadvantages of our method. In addition, the hanging drop method could not produce large-sized spheroids in 7 days, but our method was effective in producing such large spheroids. It was also found to be more efficient because increasing the number of seeded cells could reduce the number of days required for spheroid formation (data not shown). These results suggest that our method can also be applied to spheroids derived from cancer cells.

We assume that our method would be useful in cases where cells of interest are induced from pluripotent stem cells, such as iPSCs, to produce spheroids. This method involves the process of culturing in 2D for a period of time to proliferate and induce, and the process of creating spheroids to confirm the functionality of the cells of interest. Therefore, the first half of the experiment (about 7 days) was assumed to be the cell proliferation and differentiation stage and the second half was assumed to involve the production of spheroids. We assume that the greatest benefit of producing spheroids in this way is the elimination of chemical and physical stimuli. The cell line was passaged using an enzyme but was not used to produce spheroids. Since the cell line used here is passaged every 2–4 days, we minimized the effect of enzyme treatment on cells by assigning 7 days for colony formation. Thus, when increasing the seeding amount and attempting spheroidization in a short period of time, it is necessary to consider the effect of enzyme treatment.

Fluoropolymers have a high biocompatibility and inertness due to their non-polarity and strong C–F bonds. They are, thus, a safe and non-toxic material for use in the living body. Their high inertness entails advantages, such as a low adhesion in vivo, while their hampered in vivo settling is a problem for biomaterials. Their symmetrical molecular structure makes them chemically resistant and stable against most chemicals, including acids, alkalis, and organic chemicals. However, they can react slightly with molten alkali metals and their solutions as well as fluorine and chlorine trifluoride at high temperatures. Therefore, the ECM used in this study was not expected to form strong chemical bonds. In addition, the non-polar molecular structure does not allow for electrical bonds, such as hydrogen bonds. Nevertheless, the cells were able to adhere to the ECM that was transferred onto the fluoropolymer (Figure 3a,c). Hydrophobic interactions [79] and the surface charge of the fluoropolymer [78] were reported as possible mechanisms underlying this adhesion. However, we believe that van der Waals forces are instead responsible because of the weak adhesive force that allowed cells to peel off over time regardless of the ECM type. During formation of the ECM dots, the ECM floating in the culture medium settled down due to gravity and gathered on the substrate surface after the gel-like protein molecules (matrix components) were stamped. As a result, the distance between the protein molecules and the fluoropolymer surface approached zero, and subsequently, the protein molecules adhered to the fluoropolymer surface via van der Waals forces. Van der Waals forces are extremely weak compared with chemical bonds and hydrogen bonds. Therefore, when the cells proliferated and the tension between them increased, the cells began to detach from the edges of the cell colonies and were easily peeled off from the ECM-coated fluoropolymer surface upon pipetting with a weak liquid flow.

If the ECM indeed adheres via van der Waals forces, the question is why the fluoropolymer does not allow adhesion of other substances and is considered an inert material. This may be because the objects to which the fluoropolymer adheres are usually considered solids. The surfaces of solids possess microscopic irregularities, which are shown in Figure 5a. There are many gaps between the solid and the surface of the fluoropolymer. Because van der Waals forces are only generated when the distance is very short, the area ratio of the contact point to the entire interface is extremely small with solids; if the contact area/external force ratio is large, the adhesive force increases. Therefore, in many cases, it is impossible to obtain an attractive force sufficient to hold the entire solid against external forces. In contrast, when the object to be adhered is a liquid (Figure 5b), there is no gap, and the molecules are in close contact with each other. However, the molecules in the liquid move over long distances, and the van der Waals forces are not strong enough to overcome their kinetic energy and maintain the molecules at the interface. If the attractive forces are strong, the contact angle should be small, but this is not the case. Furthermore, because of the flow phenomenon, i.e., the shape of the entire liquid can freely change, it is not possible for the liquid to stay in place. In contrast, water (H2O) can remain on a glass surface because of the presence of a larger number of hydrogen bonds than in silicone oil.  The ECM used in the present study [79–84] is a gel with properties between those of a solid and liquid (Figure 5c); it can be deformed to the extent that no gaps are created between it and the fluoropolymer surface. We thus concluded that the protein molecules adhered to the fluoropolymer surface via van der Waals forces.

**Figure 5.** Application of force to the fluoropolymer depending on the state of the material placed on its surface. (**a**) Only point contact is possible due to the fine surface irregularities of solids, and thus, the area contributing to van der Waals forces is small. (**b**) Molecular motion in a spherically curled liquid as the surface tension exceeds the attractive van der Waals forces. (**c**) A gel with low surface tension and low fluidity is more densely packed, and the van der Waals forces act over a large area.

In this study, the dimensions of the spheroids were too large compared to the dimensions of the 2D matrix dots (Figure 3e). The reason is considered to be that cell lines were used here. Seven days after seeding, it was observed that the number of cell in the colony exceeded that in the 2D confluent state and proliferate in the 3D layer. If the colonies were detached by pipetting before that, the spheroid size could be controlled. In the future, the proposed method will be employed to spheroids using stem cells, and in that case, such a problem could be automatically avoided. At that time, we expect that it will also be possible to induce all cells in a planar manner at the 2D culture stage, and then suspend them while maintaining the generated cell polarity to create 3D cells using our method. In a conventional method, when a differentiation factor is added to stem cells after spheroidization, the differentiation factor may act only on the surface of the spheroid, and the induction of target cells may be insufficient in the central part. In a conventional method, when a differentiation factor is added to stem cells after spheroidization, the differentiation factors may act only on the surface of the spheroid, which can lead to insufficient induction of the target cells in the central part. In another conventional method, when a differentiation factor is added to stem cells before spheroidization, the cell substrate will be degraded during detachment via enzymatic treatment, and the cell polarity formed by 2D culture will be destroyed. Furthermore, enzymatic treatment destroys cell–cell adhesion factors, which may adversely affect the efficiency of spheroid production. Our proposed method

can solve those problems and is intended to be used for drug discovery and regenerative medicine research in the future [75,85]. We believe that this method will enable to obtain spheroids of stable quality in a simple and low-cost manner, which will contribute to the development of medicine.

#### **5. Conclusions**

We have developed a technique for generating spheroids that consists of 2D cell culture on ECM dot patterns deposited on a fluoropolymer scaffold, followed by cell colony detachment from the substrate and suspension culture. The ECM substrate was weakly bound onto the fluoropolymer surface, and thus, the adherent cell colonies were easily detached by gentle pipetting; we theorized that van der Waals forces were responsible for this weak adhesion. Spheroid formation was successful with all tested combinations of five ECM substrates (Matrigel, iMatrix, fibronectin, vitronectin, and collagen IV) and two cell lines (HepG2 and MCF7). When forming spheroids during the differentiation process of pluripotent stem cells, in contrast to conventional methods, our technique enabled spheroidization without chemical (via enzymatic treatment) or physical stimulation, and thus, a more stable spheroid production could be obtained. For spheroids such as cancer cells and cells established from tissues, however, the chemical and physical stimuli are the same as those for the conventional method. Moreover, the technique is relatively uncomplicated and is expected to be low in cost. In future studies, we will test the effect of the ECM pattern size on the spheroid size and assess differentiation induction using this scaffold to establish size controllability and a method for iPSC spheroid/organoid production. We anticipate that such spheroid/organoids can be used for drug discovery and regenerative medicine.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/app112110495/s1, Figure S1: Detailed measurement at each stage.

**Author Contributions:** Conceptualization, H.T. (Hidetaka Togo), H.T. (Hiroki Takeuchi) and M.K.; methodology, H.T. (Hidetaka Togo), H.T. (Hiroki Takeuchi) and M.K.; validation, H.T. (Hidetaka Togo), K.Y.-T. and Y.H.; formal analysis, H.T. (Hidetaka Togo); investigation, H.T. (Hidetaka Togo), K.Y.-T. and Y.H.; resources, H.T. (Hiroki Takeuchi) and M.K.; data curation, H.T. (Hidetaka Togo) and K.Y.-T.; writing—original draft preparation, H.T. (Hidetaka Togo); writing—review and editing, H.T. (Hiroki Takeuchi) and M.K.; visualization, H.T. (Hidetaka Togo); supervision, H.N.; project administration, M.K.; funding acquisition, H.T. (Hiroki Takeuchi) and M.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by JST A-STE, grant number JPMJTR20UC (to M.K.) and JSPS KAKENHI grant numbers JP18K16768, 20K18188 (to H.T. (Hiroki Takeuchi)).

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

**Informed Consent Statement:** Not applicable.

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

#### **References**

### *Article* **Application of Thermodynamics and Protein–Protein Interaction Network Topology for Discovery of Potential New Treatments for Temporal Lobe Epilepsy**

**Chang Yu 1, Edward A. Rietman 2, Hava T. Siegelmann 2, Marco Cavaglia 1 and Jack A. Tuszynski 3,4,\***

**Featured Application: This work describes the use of new methodology based on Gibbs homology analysis for the identification of potential protein targets as well as their inhibitors for the development of therapeutic options for various diseases. In the past, similar approaches have been proposed and partially validated for various types of cancer. Here, we apply the method that combines thermodynamic measures with protein–protein interaction network topology to temporal lobe epilepsy. Our results identify a number of potential therapeutic targets.**

#### **Abstract:** In this paper, we propose a bioinformatics-based method, which introduces thermodynamic measures and topological characteristics aimed to identify potential drug targets for pharmacoresistant epileptic patients. We apply the Gibbs homology analysis to the protein–protein interaction network characteristic of temporal lobe epilepsy. With the identification of key proteins involved in the disease, particularly a number of ribosomal proteins, an assessment of their inhibitors is the next logical step. The results of our work offer a direction for future development of prospective therapeutic solutions for epilepsy patients, especially those who are not responding to the current standard of care.

**Keywords:** epilepsy; systems biology; protein–protein interactions; CNS; Gibbs homology; drug targets; anti-epileptic drugs

#### **1. Introduction**

It has been estimated that about 50 million people worldwide suffer from epilepsy [1]. In 2015, about 3.4 million people had active epilepsy in the U.S. alone [2]. Epilepsy is one of the most common and most disabling neurological disorders, characterized by recurrent unprovoked excessive brain activity [3]. The current understanding of the neurophysiological mechanisms of epilepsy is largely based on extensive investigations of neuronal cells. However, glial cells have also been demonstrated to play a fundamental role in triggering seizures. Accurate diagnosis of epilepsy is challenging and therapeutic strategies span a range from single, to multiple drug courses of administration, to respective neurosurgical procedures and dietary therapy. In spite of the fact that several dozen antiepileptic drugs (AEDs) are available, there are still approximately 20–30% of patients who do not respond satisfactorily to these AEDs [3]. To develop new, effective treatments, studies focus on the development of the central nervous system (CNS) and neuronal activities in vivo to understand the causes and mechanisms of the disease initiation and progression [4]. However, to find effective treatments for patients who are refractory to current AEDs, it is

**Citation:** Yu, C.; Rietman, E.A.; Siegelmann, H.T.; Cavaglia, M.; Tuszynski, J.A. Application of Thermodynamics and Protein–Protein Interaction Network Topology for Discovery of Potential New Treatments for Temporal Lobe Epilepsy. *Appl. Sci.* **2021**, *11*, 8059. https://doi.org/10.3390/app11178059

Academic Editors: Francesco Cappello and Magdalena Gorska-Ponikowska

Received: 2 July 2021 Accepted: 27 August 2021 Published: 31 August 2021

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

**Copyright:** © 2021 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/).

necessary to study each case more precisely and individually. Hippocampal biopsy tissue of pharmaco-resistant Temporal Lobe Epileptic (TLE) patients is an extraordinarily useful substrate to study molecular mechanisms related to structural and cellular abnormalities in epilepsy. Several genes and signaling cascade alterations have already been reported in the literature based on the TLE hippocampus analysis [5]. The availability of gene expression profiling provides a detailed insight into the disease for each individual patient. Methods, such as microarray and RNA sequencing, developed in the recent decade are used to obtain genome-wide mRNA expression data. One of the advantages of using such methods is that this type of investigation is biomarker driven. Especially with whole genome-wide data, the comprehensive information on the status of each functioning unit, its interacting complex and interaction pathways is suitable for analysis, which can reveal potentially important molecular mechanisms and novel therapeutic targets. By investigating genomewide expression data, we can study the integrated results of all possible causative changes for each patient. Thus, the results of such an analysis are considered to be precise and biomarker-driven [6].

Epilepsy is a complex neuro-pathology that arises due to different etiologies, having various localizations, often occurring in conjunction with other diseases. Irrespective of what triggers the paroxysm, the electrical discharge during seizure is the common clinical manifestation for all forms of epilepsy, suggesting an underlying common molecular mechanism [7]. Although revealing the precise origin of epilepsy is still part of the ongoing investigations, it was shown that the cause may vary from de novo genetic mutations [8] to traumatic brain injury [9,10]. While gene mutations may naturally lead to altered downstream pathway behavior, brain injury was also shown to cause chronically altered gene expression signatures of genes that were linked to epilepsy [9]. Recent studies revealed some genetic causes of epilepsy, such as gene SCN1A mutations, which affect sodium channels [8] and tuberous sclerosis complex (TSC) mutations, which lead to the dysregulation of the mechanistic target of rapamycin (mTOR) pathway [11], both of which lead to epileptic conditions. Mechanisms of epileptogenesis consist of genetic and epigenetic alterations occurring in both neuronal and astroglia cells [12]. Abnormal activity of astrocytes during epileptic events has been extensively reported to play a major role due to their K+ buffering role in the extracellular milieu [13].

Clinical diagnosis of epilepsy starts with an identification of the type of seizure and proceeds with EEG and neuroimaging studies. Accurate diagnosis is pivotal for the adoption of an appropriate therapy. Anti-epileptic drugs (AEDs) represent first-line treatment for epilepsy and despite the availability of more than 20 such drugs, approximately 30% of patients do not respond to this type of therapy [12]. Among all epilepsy types, temporal lobe epilepsy (TLE) is the most common drug-resistant form of epilepsy in adults. Current AEDs mainly act by directing transmembrane ion channel function or by promoting γ-aminobutyric acid (GABA)-mediated inhibition to decrease the electrical activity of the brain [3,14]. For instance, phenytoin and lacosamide inhibit sodium channel activation; retigabine opens potassium channels; ethosuximide and lamotrigine block calcium channels; tiagabine inhibits GABA reuptake and phenobarbital and benzodiazepines enhance GABA receptors [3]. The AEDs generally focus on stabilizing and elevating the threshold of the CNS against hyperexcitability. On the other hand, the role of GABA neurotransmitters and GABA receptors in inhibiting activity in the central nervous system [15], and the role of glutamate neurotransmitters in aberrant hyperexcitability [16] have been shown to represent possible mechanisms which can be explored in order to develop new treatments. With extensive research focused on each of the above approaches, great progress has been made in deciphering epilepsy's pathophysiology at a molecular level. The better the understanding of the complex molecular mechanisms involved in the disease initiation and progression, the more unrealistic it appears that a drug affecting a single neurotransmitter receptor or an ion channel will be found to effectively treat epilepsy. Instead, combinations of pharmacological agents designed for an individual patient with a known expression profile should lead to better personalized clinical outcomes. In order to optimize such drug

combinations, sophisticated quantitative analyses of protein–protein interaction networks should be implemented involving the key biomarker proteins of interest and the results of these analysis validated experimentally. In this paper we develop such a computational modelling effort for epilepsy that is based on thermodynamic measures of protein–protein network characterization, having previously provided a general approach [17] and also applied it specifically to other diseases such as various types and stages of cancer [18–26].

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

The theoretical underpinnings for our understanding of the thermodynamics and bioenergetics of brain development started by investigating the molecular biology of human diseases from a systems and network biology perspective. These studies were developed over a several-year period [17–26]. Here, we only provide a brief summary of these approaches. The transcriptome and other -omic (e.g., proteomic, genomic, metabolomic) measures can represent the collective energetic state of a cell. By the use of the word "energetic", we mean it from a thermodynamics perspective where one uses thermodynamic functions of state such as entropy, Gibbs free energy, enthalpy, and internal energy. In particular, for open thermodynamic systems such as the human body Gibbs free energy is a suitable thermodynamic function of state, which can be computed from the so-called chemical potential for the statistical system such as a network of proteins expressed by a living cell. There, a chemical potential can be found for all interacting molecules in a cell, in particular, a chemical potential of all the proteins that interact with each other. This can be imagined to represent a rugged landscape, not dissimilar to Waddington's epigenetic landscape [27,28]. We provide mathematical expressions for the Gibbs free energy of a cellular protein–protein interaction network below.

To perform these calculations, we need input data and a method of calculating the Gibbs free energy. The method we propose uses mRNA transcriptome data or RNAseq data as a surrogate for actual measurements of protein concentration values. This assumption is largely valid since Kim et al. [29] and Wihelm et al. [30] have shown an 83% correlation between mass spectrometry proteomic information and transcriptomic information for multiple tissue types. Further, Guo et al. [31] found a Spearman correlation of 0.8 in comparing RNAseq and mRNA transcriptome from TCGA human cancer data [32]. Therefore, we have decided to use this highly correlated proxy for protein expression data in our calculations.

Given a set of transcriptome data as representative of protein concentration values, we overlay that on the graph of the human protein–protein interaction network from BioGrid [33]. This means we assign to each protein representing a node of the network, the scaled (between 0 and 1), transcriptome value (or RNAseq value). The edges in this network correspond to protein–protein interactions and they define a unique topology for a given protein–protein interaction network. As shown in our previous work [18–25], each disease studied so far is characterized by a unique network topology. From the data extracted for a given protein–protein interaction network, we compute the Gibbs free energy of each protein–protein interaction using the relation:

$$G\_i = c\_i \ln \frac{c\_i}{\sum\_j c\_j} \tag{1}$$

where *ci* is the "concentration" of the protein *i*, normalized, or rescaled, to be between 0 and 1. The sum in the denominator is taken over all protein neighbors (i.e., those that interact with it) of *i*, and including *i*. Therefore, the denominator can be considered related to a degree-entropy, although its functional form is much simplified since it does not include logarithmic terms. Carrying out this mathematical operation essentially transforms the "concentration" value assigned to each protein to a corresponding contribution to the Gibbs free energy. Thus, we replace the scalar value of transcriptome to a scalar function the Gibbs free energy. Thus, the equation represents the relationship between concentrations of proteins and the corresponding Gibbs energy.

The Gibbs free energy is a negative number, thus associated with each protein on the network that is a negative potential energy well. When plotted in 3D space where the vertical axis corresponds to the Gibbs free energy and the points in the horizontal plane represent protein coordinates, this results in a rugged energy landscape shown schematically in Figure 1. If we use what is called a topological filtration on this landscape, we essentially move a filtration plane up from the deepest energy well. As the filtration plane is moved up, larger-and-larger energetic subnetworks are captured. For convenience, we stop the filtration at energy threshold 32 meaning 32 nodes in the energetic subnetwork are retained. We call these subnetworks Gibbs-homology networks. This is not a magic number. The threshold of 32 was selected for convenience in showing networks visually. Incidentally, if we attempted to build a "network" by ranking the gene expression values, we would find disconnected nodes and not a connected network.

**Figure 1.** As the "filtration plane" moves up from the bottom, more-and-more nodes are captured in larger-and-larger energetic subnetworks shown on the right-hand side of the figure.

We now compute the Betti centrality, a topological measure, on the 32-node energetic networks as described in detail by Benzekry et al. [20]. The concept is easily explained as follows. In networks, there are holes, or rings, of various sizes. In these energetic pathways, protein–protein interaction networks, the proteins form interaction rings. In densely connected, but not fully connected, networks the rings, or holes, may consist of triangles and larger rings of interaction. To find the Betti centrality we ask ourselves the following question. Which protein when removed from the network will change the overall total number of rings the most? The total number of rings is called the Betti number. Given a network G consisting of edges, e, and vertices, v, the Betti centrality is given by

$$B(v\_i) = B(G) - B(G - \{v\_i\}).\tag{2}$$

Hence, the difference between the total Betti number *B*(*G*) and the Betti number of the network after removing node *i*, gives the Betti centrality for node labeled *i*. We compute this for all nodes in the threshold-32 energetic network. Often there will be two or more proteins in the network that have equivalent Betti centrality. We discuss this equivalence and the Betti centrality with respect to the brain region data below in this manuscript.

#### **3. Data Sources**

Patient information for TLE is available in ref. [34] using Dataset GSE63808. Note that the data from GSE63808 only include epilepsy patients and did not include healthy controls. Obviously, it would be unethical to collect temporal lobe tissue form healthy people.

Computing the Betti centrality for energy threshold 32, we find eleven proteins as the most energetically significant overall. These are ranked in the Pareto chart shown in Figure 2. (Note that Pareto ranking is a common statistical method for displaying differences in data. A Pareto diagram is a simple bar chart that ranks related measures in decreasing order of occurrence, Pareto ranking is based on the principle of non-dominated sorting also called Pareto dominance, which can clearly be seen in the present case). Since there can be one or two (sometimes three) Betti centrality nodes with equivalent energies, the number of centrality nodes in the Pareto chart adds up to greater than 131 (which is the total number of patients whose data have been accessed in this study).

**Figure 2.** Genes having the most impact on dropping Gibbs homology network complexity (threshold = 32), with counts corresponding to the number of patients having the particular gene as target.

#### **4. Results**

By using the Pareto ranking, the most important node in the network is found to be CD81. (See Figure S1 for Pareto ranking of control patients.) A Gibbs Homology network at threshold 32 in which CD81 has the highest Betti centrality is shown in Figure 3 (highlighted in yellow) and the nearest-neighbor nodes with smaller, though important Gibbs energies, are shown highlighted in green.

**Figure 3.** CD81 (yellow) and its nearest neighbors in the PPI network labeled in green.

Importantly, CD81 is a transmembrane protein that has tissue specific expression in various tissues including: tonsil, cerebral cortex, lymph node, smooth muscle, and reproductive tissues [35]. Shown in mice, CD81 regulates neuron-induced astrocytic differentiation [36]. Three alleles of the mice CD81 in seven genetic backgrounds were associated to abnormal brain development, hematopoietic system, and/or immune cells development and behavior [37]. In CD81-null mice, astrocyte and microglia cell numbers were upregulated, which lead to a 30% increase in brain size [38]. However, in human patients, deficiency of CD81 proteins was only associated to defected B cells [39]. As shown in mice and rats, expression of CD81 gene was up-regulated after seizure in three hippocampal cell layers: DGCL, CA1pyr, and CA3pyr [40,41]. In addition, the gene coding for the CD81 protein was recognized as a tumor-suppressor gene [42].

The second most frequent protein in the Pareto ranking of Betti centrality nodes is HSP90AA1 (heat shock protein 90 alpha family class A member 1). It is a stress induced isoform of the molecular chaperone Hsp90. In human, the gene has ubiquitous expression in brain and testis and 25 other tissues [43]. The expression of HSP90AA1 is known to be elevated when cellular stress is present, including in leukemia, several types of cancer, T cell under certain stimulation, and individuals suffering from chronic obstructive pulmonary disease (COPD) [44]. In rats, it was found that HSP90 protein levels in piriform cortex decreased after status epilepticus. The degradation was related to neuronal vulnerability to status epilepticus insult [45].

The third most frequent protein, excluding ribosomal proteins, according to the Pareto ranking of Betti centrality nodes is EEF1A1 (eukaryotic translation elongation factor 1, alpha 1). It plays an important role in the cellular translation process. Although EEF1A1 was not associated with neurological diseases, its isoform EEF1A2 was identified to be related. The two proteins are 92% similar in their amino acid sequences but exhibit nonoverlapping expression patterns [46]. The difference in functions that the two isoforms perform might be due to post-translational modifications [47]. In wild-type mice, expression of A2 takes over from A1 starting at 21 days after birth; at the same time, a deletion or biallelic mutation of EEF1A2 gene was found responsible for the early-onset neurological

abnormalities and early death in mice [48,49]. In human patients, EEF1A2 missense mutation was associated to early-onset epilepsy, severe intellectual disabilities and specific subtle facial dysmorphic features [50–52].