Evaluation of Cx43 gap junction inhibitors using Quantitative Structure-Activity Relationship model

1. Introduction

Gap junctions (GJs) are intercellular channels indispensable for electrical interaction between cardiac myocytes and synchronized cardiac contraction [1]. Connexin (Cx)-based GJ channels are formed of two opposing hemichannels in the contiguous cells (reviewed in [2]). Six Cx subunits compose a hemichannel with ion-selective pore. Each Cx protein has four alpha helical transmembrane domains, intracellular N- and C-termini, two extracellular loops, and a cytoplasmic loop. Cxs are expressed in all tissues except differentiated skeletal muscle, erythrocytes, and mature sperm cells, and the family of Cx genes consists of 21 members in the human genome. The prevailing connexin isoform in human cardiac tissue is Cx43 [3]. Changes in the function, expression or localization of the Cx43 are associated with a higher frequency and severity of arrhythmias and sudden death in patients with some cardiovascular diseases [3]. On the other hand, modulation of GJ function is onerous due to the shortage of specific and Cx type-selective GJ inhibitors. Most GJ inhibitors are non-specific compounds, such as antimalarial drugs, polyamines, glycyrrhetinic acid, volatile anaesthetics, arachidonic acid, cyclodextrins, anti-cancer drugs cisplatin and oxaliplatin, fatty acid amides, terpenes [4,5,6], and other drugs. In addition, most of these compounds have been shown to inhibit not only intercellular communication through GJs but also affect the membranous Na⁺, K⁺, and Ca²⁺ channels that are crucial for the generation and spread of action potential [7,8,9].

Molecular docking approaches can relatively well predict the binding affinities of compounds; however, receptor occupancy often is not directly translated to its biological effect that also depends on the intrinsic activity and intrinsic efficacy of the compound [10]. The quantitative structure-activity relationship (QSAR) is a computational modeling method used for identifying the relationship between structural characteristics of chemical compounds and their binding affinities and biological activities [11]. QSAR analysis is particularly useful in the pharmacy industry as this method helps to select substances according to their desirable biological activity, thus significantly reducing the number of substances that need to be tested in vitro and in vivo [11]. QSAR was used to develop simple, accurate, and robust models useful to identify new inhibitors [12]. For instance, in the absence of a reliable 3D receptor-bound structure more complex 3D modeling methods might face limitations [13]. Modulation of Cx43 GJ activity (activation or inhibition) is expected to be beneficial in management of the cardiac pathologies [14], such as ischemic heart disease, heart failure, hypertrophic cardiomyopathy, arrhythmogenic right ventricular cardiomyopathy, and others [15]. Therefore, the aim of this study was to evaluate known and putative Cx43 inhibitors using the molecular docking and QSAR models and compare predicted IC₅₀ (pIC₅₀) values with experimental ones (eIC₅₀) obtained by dual whole-cell patch-clamp experiments.

2. Materials and Methods

2.1. Bioinformatic analysis

The open access databases were searched for studies on the inhibitory effects of various substances on Cx43 GJs. The following substances with their half maximal inhibitory concentrations (IC₅₀) were found [14,16,17,18,19,20,21]: 2-aminoethoxydiphenyl borate (2-APB), 3,3′-methylene-bis-4-hydroxycoumarin, α-copaene, α-pinene, sabinene, carbenoxolone, digoxin, dihydrogambogic acid (DGBA), heptanol, linoleic acid, meclofenamic acid, mefloquine, quinidine, warfarin. In such a way, the dataset for Cx43 inhibition QSAR was made of 14 compounds belonging to different classes (Table 1). These classes include terpenes (α-pinene, α-copaene, sabinene), benzenes (meclofenamic acid), quinolines (mefloquine), Cinchona alkaloids (quinidine), fatty acids (linoleic acid), coumarins and derivatives (dicumarol, warfarin), triterpenoids (carbenoxolone), fatty alcohols (heptanol), cardenolide glycosides (digoxin), organoboron compound (2-APB), pyranoxanthones (DGBA). For sake of clarity, PubChem CID of all inhibitors is also provided.

Given our previous research of terpenes as Cx43 inhibitors [19] widely mentioned in literature monoterpene d-limonene [22,23,24] gained our attention as a possible Cx43 inhibitor. We decided to test this hypothesis with molecular modeling and further with dual whole-cell patch-clamp experiments.

2.2. Molecular docking

The rat Cx43 homology model was created in the Phyre2 web portal for protein modeling, prediction, and analysis [25]. Human Cx43 structure was used as a template (PDB ID 7F93) [26] for the construction of rat Cx43 homology model as a target for molecular docking. 3D molecular models for all ligands except 2-APB, digoxin, and DGBA were available in the PubChem database [27]. These 3 ligands either without clear 3D structures (digoxin and DGBA) or organoboron compound (2-APB) (no modelling parameters for B atom) were excluded from docking. Gnina version 1.0 was used for molecular docking calculations, which utilizes an ensemble of convolutional neural networks (CNNs) as a scoring function [28]. The docking mode was set to whole protein docking. All ligand output conformations except top scoring ones were excluded from further analysis automatically with program settings leaving only the most reliable one for each ligand. 3D docking image was generated with ChimeraX 1.4 software [29] and 2D docking plots were generated with LigPlot+ 2.2 software [30]. Log(pIC₅₀) from molecular docking was calculated using the equation of linear regression:

$$\log \left(pI{C}_{50}\right)=m·\mathrm{D}\mathrm{M}\mathrm{A}+n$$

(1)

where m and n are fitting coefficients and DMA is Docking Minimized Affinity.

2.3. QSAR model development

QSAR model development requires molecular descriptors calculated from their chemical structure. PaDEL-Decriptor software [31] was used to calculate 1444 molecular descriptors for those 14 compounds. Some of these descriptors had identical values for all 14 compounds like C2SP1 – amount of carbon atoms in investigated compound chemical structure containing chemical bounds to two other nearby carbon atoms and exactly one more atom of a different type had identical values for all analysed compounds so were considered useless. In total 596 molecular descriptors remained excluding those unsuitable for further modeling.

Regression analysis methods are statistical instruments widely used for the determination of relationships between molecular descriptors and biological activities of compounds [32]. Multiple linear regression (MLR) [33] for QSAR was developed by using R software for statistical computing (https://www.r-project.org/) with a Leaps package. Stepwise regression using a Forward Selection Approach [34] and also other approaches available in that package were applied. MLR was chosen to obtain relatively easily understandable QSAR models as it describes predicted biological activities in the form:

$$\log \left(pI{C}_{50}\right)=a_1{x}_1+a_2{x}_2{+}_{\dots }+a_n{x}_n+b$$

(2)

where pIC₅₀ is predicted biological activity; x_n, is molecular descriptor calculated by software; a_n is its fitting coefficient; b is the assumed log(pIC₅₀) value when all molecular descriptors are equal to 0. Given the limited size of the data set, the number of chosen descriptors for the final model (n) was chosen to be just 2 to avoid overfitting [35]. All possible models (175806 in total) with those two descriptors were created and analysed searching for one providing the most significant correlation between pIC₅₀ and eIC₅₀ values. Promising QSAR models were transferred into Microsoft Excel software for final inspection and validation. Microsoft Excel Analysis ToolPak add-in was used to compute the final values of regression analysis: R - coefficient of correlation - measured both the strength and the direction of a linear relationship; R² - coefficient of determination - provided percentage variation making it easier to compare between different models; adjusted R² helped to spot issues of overfitting. SigmaPlot 12 was also used for validation and final QSAR plot preparation.

2.4. Cell lines and culture conditions

Experiments were performed on HeLa (human cervix carcinoma, ATCC CCL-2, Manassas, VA, USA) cells stably transfected with Cx43 tagged with a green fluorescent protein (Cx43-EGFP). A stable HeLa cell line expressing Cx43-EGFP was obtained in collaboration with Dr. F. Bukauskas (Albert Einstein College of Medicine, New York, USA). The construction protocol of the vector is described elsewhere [36]. Cell line expressing Cx43-EGFP was selected using 500 µg/ml G418/geneticin (Sigma-Aldrich Co.). Cells were grown in DMEM medium containing 10% fetal bovine serum (FBS), penicillin/streptomycin mix (100 U/ml penicillin and 100 μg/ml streptomycin; Gibco Laboratories) at 37 °C and 5 % CO₂. Typically, the cells were analyzed on the second day after passage.

2.5. Electrophysiological Measurements

For electrophysiological recordings, the cells grown onto glass coverslips were transferred to an experimental chamber with constant flow-through perfusion mounted on the stage of the inverted microscope Olympus IX81 equipped with the Orca-R2 cooled digital camera. Junctional conductance g_j between contiguous cells was measured using the dual whole-cell patch-clamp technique [37]. Cell-1 and Cell-2 of a cell pair (Figure 4C) were voltage clamped independently with the patch-clamp amplifier MultiClamp 700B (Molecular Devices, Inc., USA) at the same holding potential, V₁ = V₂. By applying repeated every 10ths voltage ramps (-10 mV, 20 ms) (Figure 4D) in the Cell-1 (ΔV₁) and keeping the other constant, the junctional current was measured as the change in current in the Cell-2, I_j = ΔI₂. Thus, g_j was obtained from the ratio –I_j/ΔV₁, where ΔV₁ is equal to transjunctional voltage (V_j), and a negative sign indicates that the junctional current measured in the Cell-2 is oppositely oriented to the one measured in the Cell-1. Voltages and currents were digitized using the Digidata 1440A data acquisition system (Molecular Devices, Inc., USA) and acquired and analyzed using pClamp 10 software (Molecular Devices, Inc., USA). Patch pipettes were pulled from borosilicate glass capillary tubes with filaments. Experiments were performed at room temperature in a modified Krebs-Ringer solution (in mM): NaCl, 140; KCl, 4; CaCl₂, 2; MgCl₂, 1; glucose, 5; pyruvate, 2; HEPES, 5 (pH 7.4). Patch pipettes were filled with solution containing (in mM): KCl, 130; Na aspartate, 10; MgATP, 1; MgCl₂, 1; CaCl₂, 0.2; EGTA, 2; HEPES, 5 (pH 7.3). All chemical reagents were purchased from Sigma-Aldrich Corp. Stock solutions of d-limonene were prepared in DMSO at 100 mM concentration and later diluted with modified Krebs-Ringer solution into needed concentration (10, 30, 50 or 100 µM).

2.6. Statistical analysis

Dose–response curve obtained with different concentrations of d-limonene was fitted to a three-parameter Hill’s equation, and a concentration of the compound required to produce 50% of the maximal effect (IC₅₀) was derived using SigmaPlot 12.0 software. Data are reported as means ± SEM.

3. Results

3.1. Molecular docking of Cx43 inhibitors

Pharmacological characterization of d-limonene and other compounds was performed in silico using molecular simulation software. Molecular docking showed a common docking site for the majority of examined compounds except for sabinene, α-pinene, mefloquine, and carbenoxolone. This common docking site is a large hydrophobic furrow on the Cx43 transmembrane domains (TMs 1, 2, and 3) between neighbouring subunits (Figure 1A and Supplementary Video). Sabinene and α-pinene docked on a different site of TM3 close to extracellular loops. This hydrophobic location is suitable for docking such small highly hydrophobic compounds, but not for larger ones. Carbenoxalone docked on the Cx43 TM4 at a unique location containing two contiguous phenylalanines well fitting for its more hydrophilic part. Mefloquine docked on the pore side of TM2 reaching the N terminal domain. The exact docking locations of investigated Cx43 inhibitors are shown in 2D plots (Figure 2). From those images, it can be concluded that compounds of a highly hydrophobic nature and appropriate size are promising as Cx43 inhibitors. The docking of d-limonene is shown in 3D model (Figure 1B) which shares the same common docking site with the majority of other compounds. A fair correlation (R = 0.42) between docking minimized affinities and eIC₅₀ies was found presumably due to quite a large discrepancy between compound affinities and potencies (Figure 3A). From linear regression (Equation 1) m and n were determined equal to 0.39 and 3.99, respectively. Using these values, pIC₅₀ was calculated from docking (Table 1). Docking minimized affinity for each compound is also presented in Table 1.

3.2. QSAR modeling of Cx43 inhibitors

Using the R leaps package-based approaches the molecular descriptor pair BCUTp-1l and SM1_Dzs was selected for Cx43 inhibition QSAR modeling. BCUTp-1l is the first lowest eigenvalue of the burden matrix weighted by polarization and SM1_Dzs is the spectral moment of order 1 from the Barysz matrix weighted by I-state [38,39]. Using them the following QSAR model was developed for Cx43 inhibition:

$$\log \left(pI{C}_{50}\right)=-1.02\times BCUTp1l-0.41\times \mathrm{S}\mathrm{M}1\_\mathrm{D}\mathrm{z}\mathrm{s}+7.18$$

(3)

The developed QSAR model allowed to calculate pIC₅₀ of each compound from two calculated molecular descriptor values and compare with their eIC₅₀ies. The best QSAR model (Equation 3) achieved a very strong correlation between pIC₅₀s and eIC₅₀s: R = 0.88, R² = 0.78, and R²_adj = 0.76. Values of pIC₅₀ calculated with this model are provided in Table 1 together with literature data for comparison. The same values were used in Figure 3B where eIC₅₀ was plotted against pIC₅₀ calculated by this model. The substance with the lowest values of selected molecular descriptors (BCUTp-1l and SM1_Dzs) was heptanol and it inhibited Cx43 with the highest IC₅₀. Conversely, the compounds with the highest values of these descriptors - α-copaene and digoxin - inhibited Cx43 GJs with the lowest IC₅₀ concentration.

Figure 3. The plot of log(eIC₅₀) values versus docking minimized affinities for the docked Cx43 inhibitors (A). The plot of log(eIC₅₀) values versus log(pIC₅₀) values from the QSAR model (Equation 3) for the Cx43 inhibitors (B).

Figure 3. The plot of log(eIC₅₀) values versus docking minimized affinities for the docked Cx43 inhibitors (A). The plot of log(eIC₅₀) values versus log(pIC₅₀) values from the QSAR model (Equation 3) for the Cx43 inhibitors (B).

Table 1. Values of selected molecular descriptors (BCUTp-1l and SM1_Dzs) and values of experimental log(eIC₅₀), predicted log(pIC₅₀), and docking minimized affinities of the examined substances.

PubChem CID	Name	Log(eIC50)	BCUTp-1l	SM1_Dzs	Log(pIC50)	Docking minimized affinity	Log(pIC50) from docking
11240513	α-Pinene [19]	1.08	6.26	-0.33	0.92	-5.18	1.96
12303902	α-Copaene [19]	0.14	6.49	-0.49	0.75	-6.12	1.59
18818	Sabinene [19]	1.57	6.23	-0.64	1.08	-4.90	2.07
4037	Meclofenamic acid [18]	2.42	4.44	0.85	2.28	-6.41	1.47
4046	Mefloquine [17]	0.95	3.73	5.91	0.92	-6.66	1.38
441074	Quinidine [17]	2.60	4.91	0.92	1.77	-6.57	1.41
5280450	Linoleic acid [21]	2.00	3.72	1.19	2.88	-5.75	1.73
54676038	Dicumarol [14]	0.47	5.05	3.67	0.50	-7.54	1.03
54678486	Warfarin [14]	0.87	4.66	2.38	1.43	-6.73	1.35
636403	Carbenoxolone [18]	2.32	3.75	4.18	1.62	-7.09	1.21
8129	Heptanol [18]	3.34	3.93	0.66	2.89	-3.50	2.62
1598	2-APB [20]	1.71	4.64	0.41	2.27	-	-
2724385	Digoxin [18]	-0.86	4.26	8.45	-0.66	-	-
6857793	DGBA [16]	1.11	4.30	4.20	1.05	-	-
22311	d-Limonene	1.41	6.18	-0.50	1.07	-5.12	1.98

Table 1. Values of selected molecular descriptors (BCUTp-1l and SM1_Dzs) and values of experimental log(eIC₅₀), predicted log(pIC₅₀), and docking minimized affinities of the examined substances.

PubChem CID	Name	Log(eIC50)	BCUTp-1l	SM1_Dzs	Log(pIC50)	Docking minimized affinity	Log(pIC50) from docking
11240513	α-Pinene [19]	1.08	6.26	-0.33	0.92	-5.18	1.96
12303902	α-Copaene [19]	0.14	6.49	-0.49	0.75	-6.12	1.59
18818	Sabinene [19]	1.57	6.23	-0.64	1.08	-4.90	2.07
4037	Meclofenamic acid [18]	2.42	4.44	0.85	2.28	-6.41	1.47
4046	Mefloquine [17]	0.95	3.73	5.91	0.92	-6.66	1.38
441074	Quinidine [17]	2.60	4.91	0.92	1.77	-6.57	1.41
5280450	Linoleic acid [21]	2.00	3.72	1.19	2.88	-5.75	1.73
54676038	Dicumarol [14]	0.47	5.05	3.67	0.50	-7.54	1.03
54678486	Warfarin [14]	0.87	4.66	2.38	1.43	-6.73	1.35
636403	Carbenoxolone [18]	2.32	3.75	4.18	1.62	-7.09	1.21
8129	Heptanol [18]	3.34	3.93	0.66	2.89	-3.50	2.62
1598	2-APB [20]	1.71	4.64	0.41	2.27	-	-
2724385	Digoxin [18]	-0.86	4.26	8.45	-0.66	-	-
6857793	DGBA [16]	1.11	4.30	4.20	1.05	-	-
22311	d-Limonene	1.41	6.18	-0.50	1.07	-5.12	1.98

3.3. D-limonene dose-dependently inhibits Cx43 GJ conductance

In our earlier study, we have demonstrated that constituents of nutmeg essential oil - monoterpenes sabinene and α-pinene and sesquiterpene α-copaene - were potent and efficient Cx43 GJ inhibitors [19]. In the current study, we aimed to examine the effect of another constituent of nutmeg essential oil, suggested by QSAR modeling using the molecular descriptors, monocyclic monoterpene d-limonene, on the conductance of GJs composed of Cx43 (Figure 4 A,B) exogenously expressed in HeLa cells.

To determine the effect of d-limonene on Cx43 GJ conductance, we performed dual whole-cell patch-clamp experiments in pairs of HeLa cells expressing exogenous Cx43-EGFP (Figure 4C) applying -10 mV, 20 ms voltage ramps to Cell-1 and measuring junctional current in Cell-2 (Figure 4D). A threshold concentration of d-limonene for inhibition of g_j was 10 µM (Figure 4E). Further, applying higher concentrations (30, 50 and 100 µM) (Figure 4F–H) we found that g_j could be completely blocked by 100 µM of d-limonene. IC₅₀ value of 30 µM was derived from the fit of the experimental points to the Hill’s equation (Figure 4I). Hill’s coefficient was 2.8 suggesting more than one binding site on Cx43 GJ channel, similar to that obtained for other terpenes in our earlier study [5].

Figure 4. The effect of d-limonene on Cx43 GJ conductance. (A, B) Cx43 GJ homology model (side and top view, respectively). (C) Dual whole-cell patch-clamp measurement of Cx43-EGFP conductance in HeLa cells. (D) G_j was measured by applying repeated -10 mV, 20 ms V_j ramps which do not cause the voltage-dependent gating of Cx43 GJ channels. (E-H) Typical effects of d-limonene at indicated concentrations on Cx43 GJ conductance. (I) Dose-dependence of d-limonene effect on Cx43 GJ conductance (eIC₅₀ = 30 µM; Hill’s coefficient = 2.8).

Figure 4. The effect of d-limonene on Cx43 GJ conductance. (A, B) Cx43 GJ homology model (side and top view, respectively). (C) Dual whole-cell patch-clamp measurement of Cx43-EGFP conductance in HeLa cells. (D) G_j was measured by applying repeated -10 mV, 20 ms V_j ramps which do not cause the voltage-dependent gating of Cx43 GJ channels. (E-H) Typical effects of d-limonene at indicated concentrations on Cx43 GJ conductance. (I) Dose-dependence of d-limonene effect on Cx43 GJ conductance (eIC₅₀ = 30 µM; Hill’s coefficient = 2.8).

4. Discussion

Cardiac remodeling, which involves structural and electrical changes in the heart, may be impacted by altered expression and localization of Cx43 GJs. Decreased expression of Cx43 proteins and heterogeneous arrangement of channels can impair cardiac conduction and lead to supraventricular or ventricular arrhythmias [40]. Cx channels are promising pharmacological targets because inhibitors of Cx channels could be useful for treating not only arrhythmias, but also other communication-dependent diseases affecting other body systems. The importance of Cx43 has been well established, particularly in the heart where a lack of Cx43 leads to abnormal cardiac development and death at birth [41]. Abnormalities in Cx43 organization and regulation have also been linked to myocardial ischemia [42]. Preservation of Cx43 in the intercalated discs joining cardiomyocytes may not only support conduction velocity but also impact the total normoxic Cx43 interactome [43]. So, Cx43 is a considerable drug target, especially during heart ischemia and reperfusion [14,44,45]. For example, Cx43 inhibitor digoxin (analyzed in QSAR here) isolated from Digitalis lanata is well known in cardiology [46]. It is used to treat both irregular heartbeat [47] and heart failure [48], but its side effects like gynecomastia also are significant which could be explained through its chemical similarity to estrogen [49]. Another Cx43 inhibitor quinidine also is a popular antiarrhythmic drug [50]. Another examined compound carbenoxolone is used for the treatment of peptic, esophageal and oral ulceration and inflammation; however, it also has been shown in humans to slow the myocardial conduction [51]. An anti-malarial drug mefloquine may lead to complete heart arrest [52] and number of neuropsychiatric effects including suicide [53]. More recently, Cx inhibiting peptides (antiarrhythmic peptides AAP10, ZP123; GAP-134; RXP-E; and the Cx43 mimetic peptides Gap 26 and Gap 27) were suggested for treatment of arrhythmias in patients with ischaemic cardiomyopathy [54]. On the other hand, peptides underperformed as drug candidates due to unfavorable characteristics mainly regarding their pharmacokinetic behavior, including plasma stability, membrane permeability, and circulation half-life [55]. The discovery of new modulators of GJ channel function lacking similarity to steroid hormones to avoid side effects is of interest for human health [14,44,45]. Also, a long-standing challenge in the study of GJs is the lack of specific, high-affinity activators and inhibitors of GJ channels [56,57]. Therefore, it is important to predict in silico which substances could effectively modulate GJ conductance. This also could serve as an innovative approach in repurposing of licensed drugs with predicted new Cx inhibitory properties for other communication-dependent illnesses.

The results obtained in molecular docking considering minimized affinity had much worse correlation with log(eIC₅₀) (R was 0.42) when compared to QSAR modeling (R was 0.88) so it might be concluded that QSAR outperforms docking. Fair correlation in case of molecular docking can be explained by often happening discrepancy between affinity and potency of the compounds [58].

All top 3 Cx43 inhibitors docked into the same common docking site suggesting that this site should be explored further for more potent inhibitors. Thus, even if molecular docking accuracy is low it can be used in combination with QSAR. Using the developed QSAR model, the limonene log(pIC₅₀) was calculated equal to 1.07 (Equation 3), and an experimental value of 1.41 was later determined by our patch-clamp experiments. Therefore, this new Cx43 inhibitor could be added the current GJ inhibitor nomenclature.