Model-driven engineering of supramolecular buffering by multivalency.

A supramolecular system in which the concentration of a molecule is buffered over several orders of magnitude is presented. Molecular buffering is achieved as a result of competition in a ring-chain equilibrium of multivalent ureidopyrimidinone monomers and a monovalent naphthyridine molecule which acts as an end-capper. While we previously only considered divalent ureidopyrimidinone monomers we now present a model-driven engineering approach to improve molecular buffering using multivalent ring-chain systems. Our theoretical models reveal an odd-even effect where even-valent molecules show superior buffering capabilities. Furthermore, we predict that supramolecular buffering can be significantly improved using a tetravalent instead of a divalent molecule, since the tetravalent molecule can form two intramolecular rings with different "stabilities" due to statistical effects. Our model predictions are validated against experimental 1H NMR data, demonstrating that model-driven engineering has considerable potential in supramolecular chemistry.

The high level of complexity found in biochemical systems often necessitates the synergy of a combined experimental and theoretical study (1, 2). Moreover, a well-established approach in the fields of systems and synthetic biology is to develop novel functionalities by modeling the required molecular mechanisms before any experimental work is performed, i.e., model-driven engineering (3). Thus, model-driven engineering is defined here as the intensive utilization of computing and informatics with the aim to assess the viability of novel molecular designs and to extract design rules, thereby decreasing the time required for experimental work. The development of synthetic supramolecular systems is currently at a level of complexity that requires a similar synergistic experimental and theoretical treatment (4–6). Analogously, theoretical descriptions are beginning to approach the required level of predictive accuracy required for model-driven engineering (7).

Multivalency is a ubiquitous phenomenon in biochemical systems that is associated with high binding affinity, increased selectivity, as well as ultrasensitivity (8–12). Since these properties can be of invaluable use in unprecedented molecular engineering approaches, multivalency is often applied in synthetic systems (13, 14). For example, research has shown that multivalent medication can have much lower toxicity while simultaneously having higher medical efficacy (10, 15). More recently, multivalency has been recognized as a key molecular driving force in the formation of membraneless organelles in living cells (16). These phase-separated cellular bodies are organized by dynamic multivalent interactions between proteins and RNA scaffolds and offer a compartmentalized liquid environment that promotes specific enzymatic reactions due to high local concentrations and insulates these reactions from competing substrates (17).

Theoretical and experimental studies of multivalent systems have revealed several design parameters that are critical in obtaining effective multivalent constructs. Next to the binding affinity, linker flexibility plays an important role: Rigid linkers require extremely precise ligand positioning to obtain high binding affinities and selectivity, while flexible linkers offer more freedom in molecular design at the cost of lower affinity and selectivity (18–20). Furthermore, additional competing equilibria can be used to enhance binding selectivity or to steer an assembly toward a preferred state (5, 21).

Most, if not all, multivalent biochemical systems are based on dimerization by specific host–guest binding with minimal host–host and/or guest–guest interactions. As a result, most studies in multivalency focus on those types of systems. However, with the rise of supramolecular chemistry, an increasing number of multivalent constructs that have self-associating groups are becoming available. The addition of self-association broadens the behavior of the multivalent constructs to include more possibilities for intramolecular cyclization. In divalent homodimerizing systems, this can have various interesting effects such as mechanically induced gelation, entropy-driven polymerization, or light-switchable gelation (22–24). Self-associating constructs with higher valencies are reported less often and are typically used for their gelation properties, where cyclization leads to less “effective” gelation, or as polymer glasses (25–28).

Recently, our group reported on a two-component supramolecular buffering system based on a self-associating divalent ureidopyrimidinone (UPy) molecule and a monovalent naphthyridine (NaPy) molecule that undergoes dimerization with UPy chain ends (29, 30) (Fig. 1). The term “supramolecular buffering” refers to the insensitivity of the concentration of a free component with respect to changes in the total concentration. In this case, the buffered free component is the monovalent NaPy which acts as a stopper molecule for the divalent UPy. The buffering of NaPy originates from the competition between cyclization of the divalent UPy and end-capping of linear oligomers by NaPy. When the total concentration of divalent and monovalent components is changed simultaneously, this leads to a buffering plateau in which the free NaPy concentration is constant over several orders of magnitude (Fig. 1). Since the effectiveness of the buffering is controlled to a large degree by the cyclization tendency, we hypothesized that multivalent constructs with higher valencies might lead to improved buffering. Therefore, we present a systematic study in which we investigate how multivalency affects supramolecular buffering using a model-driven engineering approach.

Fig. 1.

(A) Equilibrium between a divalent UPy and monovalent NaPy, which governs supramolecular buffering. UPy–UPy and UPy–NaPy dimerization is possible due to fourfold hydrogen bonding after tautomerization of the UPy moiety. (B) Supramolecular buffering of the free NaPy concentration (CNaPy free) versus the total concentration. Buffering occurs in region II due to competition between cyclization of the divalent UPy and end-capping of linear species by NaPy. Adapted with permission from ref. 29. Copyright 2015 American Chemical Society.

Results

To analyze how multivalency affects supramolecular buffering, we expanded our previous model describing supramolecular buffering by divalent UPy monomers (29). To this end, we developed models that describe ring–chain competition of tri- and tetravalent UPy monomers, followed by the inclusion of NaPy dimerization with UPy chain ends. Models for the dimerization of monovalent molecules have already been established (31). While there is a multitude of theoretical models available that describe ring–chain equilibria of tri- or tetravalent molecules with host–guest binding groups, self-association has, to the best of our knowledge, not been included (10, 15, 32–34). Furthermore, even for noncyclizing tri- and tetravalent molecules, it is analytically intractable to include aggregates with high degrees of polymerization (DPs) due to the exponential increase in the number of molecular species as a function of DP. Moreover, the inclusion of cyclization further increases the number of distinct species. However, since the stopper molecule will limit the formation of larger assemblies, the expected DP will remain low at intermediate concentrations as studied here. Thus, mass balances for the tri- and tetravalent molecules were constructed up to a DP of four for the multivalent molecule ().

The Jacobson–Stockmayer theory, which describes the polymerization and cyclization of divalent molecules in reversible covalent polymerizations, forms the basis for the constructed models (35). The theory has been refined by allowing finite intermolecular binding constants (K), i.e., supramolecular contacts instead of covalent bonds (36). In Jacobson–Stockmayer theory, cyclization is taken into account via the effective molarity (EMi), which is the experimentally measured cyclization tendency of a chain consisting of i divalent molecules. When the value of EMi for any i is known, the remaining values can be predicted by assuming that the linker follows Gaussian chain statistics, i.e., the linker is strainless. If the linker is not strainless, which is sometimes observed experimentally for relatively short linkers (<30 atoms), the behavior can be described using multiple EMi values (37).

In the theoretical models for the tri- and tetravalent molecules, cycle stabilities are calculated by assuming that the ring–closure equilibrium constant of a chain with i tri- or tetravalent molecules is equal to that of a chain of i divalent molecules, while differences in statistical factors are taken into account (Fig. 2 ; for full model details, see ). Statistical factors are used to enumerate the number of ways a reaction can proceed forward and backward, thereby accounting for the degeneracy of assembled states (32). The subsequent inclusion of NaPy dimerization with UPy chain ends leads to a considerable increase in the number of species as a result of combinatorial complexity (38) (). The input parameters for the models are the intermolecular UPy–UPy and UPy–NaPy equilibrium binding constants (KUPy–UPy and KUPy–NaPy), the effective molarity of the monomeric ring (EM1), and the ratio of NaPy to multivalent UPy (f).

Fig. 2.

(A and B) Selection of species used in the tri- and tetravalent UPy models. Partially cyclized species and species larger than dimeric are omitted for clarity. Statistical factors are based on the reference reactions of inter- and intramolecular UPy–UPy dimerization shown in the squares. (C) Simulated stopper titration showing the predicted free NaPy concentration versus the total NaPy concentration (solid lines) and the maximum twofold sensitivity (Inset). (D) Simulated equimolar dilution showing the predicted free NaPy concentration versus the total multivalent UPy and NaPy concentration (main) and the fraction of cycles during the dilution (Inset). The dashed lines indicate unreliable model predictions, where the fraction of species that have the highest degree of polymerization included in the model is larger than 30%. The input parameters for the predictions in C and D are CmultiUPy = 10 mM, KUPy–UPy = 6 × 107 M−1, KUPy–NaPy = 5 × 106 M−1, and EM1,divalent = 10 mM.

Using estimated values for the parameters that are appropriate for the multivalent UPy and NaPy system in chloroform, both stopper titrations and 1:1 dilutions were simulated (Fig. 2 ). In the stopper titration simulations, the free NaPy concentration is calculated as a function of the total concentration of NaPy while the concentration of multivalent UPy is constant. The multivalent UPy concentration is set equal to the EM of the divalent molecule to ensure that no linear species are present before NaPy is added. The simulated titration curves, displaying the free concentration of NaPy as a function of the total NaPy concentration, reveal a shallow slope at low NaPy concentrations and a steep transition at the equivalence point for all valencies (Fig. 2). Interestingly, at low NaPy concentrations both the mono- and trivalent UPy titration curves have a slope of unity, while the di- and tetravalent UPy curves have lower slopes, the latter indicative of buffering of free NaPy by cyclic species (29). Furthermore, the slope of the tetravalent UPy titration curve at low NaPy concentrations is lower than that of the divalent. This lower slope is attributed to the fact that the two intramolecular cycles of the tetravalent UPy monomer have different stabilities. The difference in stability is not due to any cooperativity, but solely due to a difference in statistical factors for ring formation (). Thus, during the stopper titration, the binding of NaPy stopper to the cyclized tetravalent molecule proceeds in a two-step process, opening the cycles one by one (). This effect suggests that a multivalent construct with an even higher valency might lower the slope to zero, assuming it would be able to dissolve. As expected, the simulated trivalent UPy titration curve overlaps with that of the monovalent UPy at low concentrations, since the first UPy–UPy contact to break upon addition of NaPy is the intermolecular contact (5). Thus, at low concentrations, the trivalent UPy acts as a monovalent UPy dimer with two additional cycles.

At NaPy concentrations around the equivalence point, all simulated curves show a sharp increase in the free NaPy concentration, similar to the ultrasensitive response observed in molecular titration (11) (Fig. 2). The magnitude of the response in this concentration regime can be characterized by the twofold sensitivity, which is defined as the change in output (free NaPy concentration) when a twofold change in input is applied (total NaPy concentration) (11). Interestingly, the titration curves corresponding to the mono- and divalent UPy constructs have similar twofold sensitivities, while the titration curves of the tri- and tetravalent constructs have increasing sensitivities, indicating sharper transitions (Fig. 2, Inset). The twofold sensitivities are also dependent on the concentration of the multivalent UPy constructs, as it determines the fraction of cyclic species at the start of the titration and the resulting degree of competition between cyclization and end-capping (). Thus, multivalency can be used to generate sharper transitions and to improve supramolecular buffering in this system.

For the 1:1 dilution simulations, both the concentrations of multivalent UPy and NaPy are changed simultaneously, keeping the ratio between the two constant at a value of unity (Fig. 2). Previously, we have shown that for a divalent UPy construct, the concentration of free NaPy is independent of the total concentration of NaPy and UPy over a broad concentration regime (29). The simulations reveal that both the tri- and tetravalent UPy constructs should display a broad plateau indicative of supramolecular buffering (Fig. 2). Interestingly, while the supramolecular buffering of the trivalent UPy construct only occurs at higher concentrations (compared with the divalent construct), the tetravalent UPy construct buffers at the same concentrations as the divalent construct but yields a broader buffering plateau. The inferior buffering by the trivalent construct is not entirely unexpected, since competition between cyclization and end-capping is a key requirement for supramolecular buffering. Because the first association of the NaPy stopper with the trivalent molecule will initially disrupt any intermolecular UPy–UPy contacts, buffering is expected to occur only at higher concentrations when cycles are opened (Fig. 2, Inset). The superior buffering by the tetravalent molecule is attributed to the fact that it lacks a critical concentration above which only chains are present, as is the case for divalent ring–chain equilibria (29, 35, 36). Instead, the formation of intermolecular contacts does not prevent cyclization of the remaining UPy moieties, allowing further competition between cyclization and end-capping and a continuation of the buffering plateau (Fig. 2, Inset). Therefore, it is expected that multivalent constructs with even higher valencies will show a similarly extended supramolecular buffering plateau.

To validate the model predictions, a library of multivalent UPy molecules was synthesized with valencies ranging from mono- to tetravalent (Fig. 3). In an effort to exclude influence from steric repulsion on the linker flexibility and subsequently the cyclization tendency of multivalent UPys 2–4, the di- and trivalent UPys (2 and 3) were equipped with methyl groups at the central branching position. While this approach does not completely exclude variations in the linker flexibility, it does prevent any additional attractive interactions between linker segments by avoiding heteroatoms.

Fig. 3.

(A) Molecular structures of multivalent UPys 1–4, NaPy 5, and DPU 6, used in this study. (B) 1H NMR data (squares) and model predictions (lines) of the free NaPy concentration in dilution experiments using mixtures containing divalent UPy 2 and f equivalents of NaPy 5. (C and D) 1H NMR data (squares) and model predictions (lines) of the free NaPy concentration in dilution experiments using mixtures containing f equivalents of NaPy 5 and one equivalent of either trivalent UPy 3 (C) or tetravalent UPy 4 (D). Error bars denote 95% confidence intervals based on assumed relative SDs of the mass, volume, and NMR integral (1, 1, and 5%, respectively) (29). All 1H NMR experiments were performed in CDCl3.

The synthesis of the multivalent UPy library was initially focused on the tetravalent UPy molecule, since its solubility was expected to present an obstacle. Several attempts were made to synthesize tetravalent UPys with urethane groups in the linker and with varying substituents on the 6 position of the pyrimidinone ring (methyl, ethylpentyl, adamantyl). However, those approaches all resulted in precipitate formation during the final reaction step, which could not be redissolved. Thus, to maximize the solubility, the urethane in the linker was replaced with an ester group and UPy groups with an ethylester on the 5 position were employed. Gratifyingly, this successfully yielded tetravalent UPy 4 which has a reasonable solubility (∼4 mM in CHCl3). The synthesis of multivalent UPys 1–3 was performed in a similar manner as that of tetravalent UPy 4 and proceeded without further problems ().

To determine the correct model parameters to be used in the validation of the tri- and tetravalent UPy models, KUPy–NaPy and EM1 were determined experimentally. The value of KUPy–UPy is fixed at the reported value for 6-methylureidopyrimidinone groups in CHCl3 (KUPy–UPy = 6 × 107 M−1) (29, 39). The correlated value of KUPy–NaPy was determined by first measuring 1H NMR spectra of equimolar mixtures containing monovalent UPy 1 and NaPy 5, followed by fitting of the measured distribution of UPy–UPy and UPy–NaPy contacts with a simple binding model that includes self-association of the UPy groups and UPy–NaPy dimerization [KUPy–NaPy = (3.1 ± 0.2) × 106 M−1 in CHCl3; See ]. The determination of the EM of the divalent molecule was performed by measuring concentration-dependent 1H NMR spectra of divalent UPy 2 and subsequently fitting the data with a ring–chain model for divalent molecules (EM1 = 5.3 ± 0.3 mM; see ) (36).

Gratifyingly, using the optimized values of KUPy–NaPy and EM1 to predict the buffering in a 1:1 dilution experiment of NaPy 5 and divalent UPy 2 mixtures results in an excellent prediction of the free NaPy concentration over a broad concentration range (Fig. 3). Thus, with the parameters for dimerization and cyclization determined, the model predictions of the tri- and tetravalent models were tested. Various dilution experiments on mixtures of NaPy 5 and either trivalent UPy 3 or tetravalent UPy 4 were performed using 1H NMR spectroscopy while keeping f, the ratio between NaPy and multivalent UPy, constant (Fig. 3 ). While the model predictions completely overlap with the measured free NaPy concentration within the confidence bounds, we note that at low fractions of free NaPy the measurement of free NaPy concentration is not reliable. The free NaPy concentration is calculated indirectly from the total NaPy concentration and the concentration of UPy–NaPy contacts, which leads to an increase in uncertainty at low free NaPy fractions (29). The concentration of UPy–NaPy contacts is obtained from the UPy N-H resonances in the 1H NMR spectra. To better validate both models, two separate global fits of the UPy N-H resonances of 3 and 4 were performed with three free parameters (KUPy–UPy, KUPy–NaPy, and EM1; See ). Since the UPy speciation results directly from the 1H NMR spectra it is more reliable than the indirect calculation of the free NaPy concentration (). Gratifyingly, the values of the optimized parameters were close to the experimentally determined parameter values.

To further validate the tetravalent UPy model, several titrations with N,N’-di-2-pyridylurea (DPU) 6 were performed on equimolar mixtures of tetravalent UPy 4 and NaPy 5 in CDCl3 (Fig. 4 and ). DPU 6 selectively binds to NaPy 5 due to its complementary ADDA hydrogen bonding array, effectively sequestering NaPy from the mixture (40). Furthermore, DPU 6 has no interaction with UPy groups, since the UPy groups cannot tautomerize to the complementary DAAD configuration. 1H NMR spectra obtained during the titration showed that the UPy–NaPy contacts were disrupted, which is in line with NaPy sequestration by DPU, and that the fraction of UPy–UPy contacts in monomeric rings increased, which is consistent with the dilution of the tetravalent UPy (Fig. 4 ).

Fig. 4.

(A) Representation of the titration of equimolar mixtures of tetravalent UPy 4 and NaPy 5, using DPU 6 as the titrant. The addition of DPU causes UPy–NaPy contacts to break in favor of UPy–UPy and NaPy–DPU contacts. The titrant solution contains only DPU, which causes dilution of tetravalent UPy 4 and NaPy 5 during the titration. (B) Fit (lines) of 1H NMR data (markers) of the DPU titration (CtetravalentUPy,0 = CNaPy,0 = 1 mM in CDCl3). (C) UPy speciation. (D) NaPy contacts during the titration. Only NaPy–UPy contacts could be deduced from the 1H NMR spectrum. (E) NaPy speciation. In the speciation plots (C and E), species with a DP >3 of the tetravalent UPy are omitted for clarity. The colors of the lines correspond to the DP of the tetravalent UPy (1: blue, 2: red). Error bars denote 95% confidence intervals based on assumed relative SDs of the mass, volume, and NMR integral (1, 1, and 5%, respectively).

To validate the titration data against model predictions, the tetravalent UPy model was adapted to include NaPy–DPU dimerization and DPU self-association (). Since the reported binding constants for DPU were only approximately determined, a global fit of the titration data was performed using two free parameters (KDPU and KDPU-DPU). The values for KUPy–NaPy and EM were set to those determined by the reference experiments with monovalent UPy 1, NaPy 5, and divalent UPy 2, vide supra. The titration data could be fitted well, and the optimized parameters are in close correspondence with the reported approximate values (Fig. 4 and ).

Discussion

Both the DPU titration and the data of multivalent UPy–NaPy mixtures show that the models for the tri- and tetravalent molecules can sufficiently describe the behavior of the multivalent molecules in the presence of NaPy stopper. Therefore, the models are correct in predicting the inferior buffering of the trivalent UPy and the superior buffering of the tetravalent UPy molecule, compared with the divalent construct. The good agreement between model predictions and experimental results show that model-driven engineering is an outstanding strategy to investigate new molecular topologies. While this improved supramolecular buffering system does not yet approach the performance of pH buffers in titration experiments, where slopes of zero can be obtained, we do show that multivalency can be used to improve both the capacity of the buffer and the sensitivity of the response.

Our study on the effects of multivalency on supramolecular buffering revealed an odd–even effect, where the buffering by molecules with odd-numbered valencies is significantly inferior to molecules with even valencies. Curiously, a similar odd–even effect was found in catalytic activity of multivalent dendrimers equipped with catalysts capable of both single and double-site catalysis (41). Supramolecular buffering can be substantially improved by employing a tetravalent molecule, as it is able to form two intramolecular rings with different stabilities due to statistical factors. Furthermore, we show that multivalency, while mostly employed to generate sharper responses, can also generate systems that are insensitive to changes in concentration.

The present system might be further developed by utilizing the tetravalent molecules described here analogously to multivalent protein and RNA constructs that form phase-separated cellular bodies (17, 42). Cellular bodies use phase separation to buffer components, isolate incompatible substrates or catalysts, and promote specific reaction rates by changes in local concentrations. Therefore, such an approach may provide a next step in the construction of artificial cells while simultaneously providing a fundamental framework for the effects of phase separation.

The multivalent constructs could also be incorporated in chemical reaction networks, as combining multivalency and catalysis could lead to increased control over the reaction rate. It would allow for sharper switching between the on and off state and higher rates of catalysis due to increased local concentration. The kinetics of multivalent catalysts and multivalent substrates have been investigated in detail (41, 43), and their incorporation in chemical reaction networks could shed more light on analogous molecular mechanisms in biochemical pathways.

Materials and Methods

Simulations were performed using the MATLAB software package (R2016a, Version 9.0.0341360; MathWorks) along with its optimization, curve fitting, and symbolic math toolboxes. Where appropriate, mass balances were analytically solved using the Mathematica software package (Version 9.0.1.0; Wolfram Research, Inc.). Otherwise, mass balances were solved numerically using either the fzero or fsolve function included in MATLAB. Nonlinear least-squares optimizations were performed using the lsqcurvefit function from MATLAB’s optimization toolbox. This function uses the Levenberg–Marquardt method to minimize the residual sum of squares. For each optimization 1,000 fits were performed. Initial parameters for the fits were distributed using latin hypercube sampling (implemented in the lhsdesign function), which ensures a uniform distribution in multidimensional parameter space so that the global optimum can be obtained. The optimization with the lowest squared two-norm is used as the best fit, while optimizations with a squared two-norm within 5% of the best fit are considered equally good fits.

All solvents, except deuterated solvents, were obtained from Biosolve; all other chemicals were purchased from Sigma-Aldrich unless noted otherwise. Deuterated chloroform was obtained from Cambridge Isotope Laboratories. Dry CHCl3 was obtained by adding oven-dried molecular sieves (4 Å) at least 48 h before the measurements.

NMR spectra were taken on a 400- or 500-MHz Varian spectrometer and results were processed using Mestrenova software. Proton chemical shifts are reported in ppm downfield from tetramethylsilane (TMS) and carbon chemical shifts in ppm downfield of TMS using the resonance of the deuterated solvent as internal standard. Abbreviations used are b, broad; d, doublet; d-d, double doublet; p, pentet; q, quartet; s, singlet; and t, triplet. Flash column chromatography was performed on a Biotage Isolera Spektra One Flash Chromatography system using KP-Sil Silica Gel SNAP columns. Mass spectrometric characterization was performed on a Bruker Autoflex Speed MALDI-TOF spectrometer. Sample preparation was performed using 2-[(2E)-3-(4-tert-Butylphenyl)-2-methylprop-2-enylidene]-malononitrile or α-cyanohydroxycinnamic acid as the matrix. Recycle gel permeation chromatography (GPC) was performed using a Shimadzu system equipped with a Jai-Gel 2.5 H and a Jai-Gel 2 H column in series employing UV-Vis detection at 275 nm and 325 nm.

Synthetic procedures, supporting simulations, parameter determinations, models for the tri- and tetravalent UPys, and the fits for model validation are detailed in .