Self-Supporting Hyaluronic Acid-Functionalized G-Quadruplex-Based Perfusable Multicomponent Hydrogels Embedded in Photo-Cross-Linkable Matrices for Bioapplications

Dynamic G-quadruplex supramolecular hydrogels have aroused great interest in a broad range of bioapplications. However, neither the development of native extracellular matrix (ECM)-derived natural biopolymer-functionalized G-quadruplex hydrogels nor their use to create perfusable self-supporting hydrogels has been explored to date, despite their intrinsic potential as carrier vehicles of therapeutic agents, or even living cells in advanced regenerative therapies, or as platforms to enable the diffusion of nutrients and oxygen to sustain long-term cell survival. Herein, we developed a dynamic co-assembling multicomponent system that integrates guanosine (G), 3-aminophenylboronic acid functionalized hyaluronic acid (HA-PBA), and potassium chloride to bioengineer strong, homogeneous, and transparent HA-functionalized G-quadruplex hydrogels with injectable, thermo-reversible, conductive, and self-healing properties. The supramolecular polymeric hydrogels were developed by hydrogen bonding and π–π stacking interactions between G coupled via dynamic covalent boronate ester bonds to HA-PBA and stabilized by K+ ions, as demonstrated by a combination of experiments and molecular dynamics simulations. The intrinsic instability of the self-assembled G-quadruplex structures was used to bioengineer self-supporting perfusable multicomponent hydrogels with interconnected size and shape-tunable hollow microchannels when embedded in 3D methacrylated gelatin supporting matrices. The microchannel-embedded 3D constructs have shown enhanced cell viability when compared to the bulk hydrogels, holding great promise for being use as artificial vessels for enabling the diffusion of nutrients and oxygen essential for cell survival. The proposed approach opens new avenues on the use of ECM-derived natural biopolymer-functionalized dynamic G-quadruplex hydrogels to design next-generation smart systems for being used in tissue regeneration, drug screening, or organ-on-a-chip.


Extended Discussion
Following the typical structural assessment of DNA-based G-quadruplexes, the impact of the different cation binders was evaluated with resort to the Root-Mean-Square-Deviations (RMSD) of the G4's, plotted in Figure S3 throughout the MD simulation time. The average values for the concatenated final 10 ns of each MD run are given in  Table S2 and Figure S4), due to the bulky HA backbone, yields larger values than in DNA-based G-quadruplexes (ca. 17 vs 11 Å). 1 However, the assessment of the G4's ROG leads to values close to 6 Å, being also insufficient to discriminate between the different cations as G-quadruplex binders.        3.58 ± 0.27 [3.15 : 5.02] a) N = 90000 frames (10000 frames × 3 MD runs × 3 intervals); and b) When considering the best-preserved pair of G4´s the average distance is 3.86 ± 0.56 Å, in the 3.12 to 6.95 Å range. Scheme S1. Putative hydrogen bonding interactions between the guanine moieties of A, with the corresponding atom numbering scheme.

Computational Methods
The The force field parametrisation of the boron bonding terms was carried out using a smaller fragment of A, Afrag, sketched in Scheme S3.
Scheme S3. Afragfragment of A used to derive the force field parameters.

Quantum Calculations
The MD simulations were preceded by the development of force field parameters comprising bonding terms for boron and derivatisation of restrained electrostatic potential (RESP) atomic charges, which were performed with Gaussian 16, 8 as detailed below.

Classical force field calculations
The MD simulations were carried out with Amber20, 9 using parameters taken from the second generation of the Generalized Amber Force Field (GAFF2), 10,11 except those involving the boron centre, coupled with RESP atomic charges. 12 The water molecules were described with the TIP3P model, 13 the monoatomic cations were described with a discrete charge of +1 or +2 and van der Waals parameters developed for the TIP3P water model, 14,15 whilst NH4 + was described with GAFF2 and standard RESP charges. [10][11][12] The post-processing of trajectory files to obtain the structural data was performed with cpptraj. 16

Boron force field parameters
The boron bonding terms were derived using the VFFDT software 17 and the paramfit utility, 18 available within the Amber20 software package. This technical approach has been successfully used in the parametrisation of boronate derivatives. 19,20 The boron van der Waals parameters were directly taken from reference. 21 The bond length and bond angle terms were estimated for Afrag using the Seminario method, as implemented in VFFDT, 17  The dihedral angle terms centred at the oh-B and ca-B bonds were calculated with paramfit, using the following workflow. Afrag underwent two independent Energy Surface scans through DFT single point energy calculations, with the dihedral angles centred at the oh-B or ca-B bonds being increased by 2° 180 times, resulting in a 360° scan around each bond.
Subsequently, the dihedral force field parameters were obtained using the genetic algorithm of paramfit in which the MM energies were fitted to the quantum data. The dihedral parameters around the os-B bonds were adapted from the X-os-c3-X parameters available in GAFF2.
Unless otherwise specified, all remaining dihedral bond angle parameters as well as the improper terms were adapted from GAFF2, considering that the boron is equivalent to a tetrahedral c3 carbon atom type, as summarised in Table S6.

RESP atomic charges of A
A fragment-based approach was employed for the calculation of the RESP atomic charges,

General MD simulation methods
An initial G4 of four A molecules surrounding one K + cation was built and minimised in the gas phase until the convergence criterion of 0.0001 kcal mol -1 was achieved. Afterwards, the structure of the minimised G4 was repeated three times over itself, affording a G-quadruplex assembled by four K + coordinated to the guanines' oxygen atoms. This model structure was minimised in the gas phase by MM and used as a starting binding scenario for the subsequent MD simulations, regardless of the cation. Afterwards, this structure was solvated in a truncated octahedron box with 15400 TIP3P model water molecules. For each system, a variable number of cations was added to neutralise its net charge.
Each solvated system was equilibrated under periodic boundary conditions using the following multistage protocol. The system was relaxed by MM minimization of solvent molecules and by keeping all solutes fixed with a positional restraint of 500 kcal mol -1 Å -2 . The restraint was then removed, and the entire system was allowed to relax. These two minimization stages comprised an initial set of 10000 steepest descent algorithm steps, followed by 10000 steps of conjugated gradient algorithm. The equilibration proceeded with heating up the system to 300 K for 100 ps using an NVT ensemble and a weak positional restraint (10 kcal mol -1 Å -2 ) on the G-quadruplexes and coordinated cations. Afterwards, maintaining the positional restraints, each system's density was allowed to equilibrate in an NPT ensemble at 1 atm for 1 ns, at the same temperature, followed by an NPT data collection run of 50 ns without any positional restraints. The collection run's trajectory frames were saved every 1 ps. Three independent runs were performed for each system. The CUDA version of the PMEMD executable was used for the simulation of all solvated systems. 24,25 The bond lengths involving all bonds to hydrogen atoms were constrained with the SHAKE algorithm, allowing the usage of a 2 fs time step. 26 The Particle Mesh Ewald (PME) method was used to treat the long-range electrostatic interactions. 27 The non-bonded van der Waals interactions were truncated with a 10 Å cut-off. Figure S6. 1 H NMR spectra of G-quadruplex hydrogel at pH 4, 7.4, and 10 in D2O, 300.13 MHz; Ha and Hb refer to reacted and unreacted G, respectively.