The Stupp laboratory has developed a class of peptide amphiphiles that spontaneously form high-aspect-ratio, filamentous nanostructures in water (Hartgerink et al., 2001). These peptide amphiphile molecules typically consist of a hydrophobic hydrocarbon tail conjugated to a sequence of amino acids with the propensity to form β-sheet hydrogen bonds and charged amino acids for solubility. The headgroups of these peptides typically contain amino acid residues with weakly acidic or weakly basic side chains (Hendricks et al., 2017). The competition between the hydrophobic interactions of the alkyl chain, the packing of the alkyl chain and headgroup, electrostatic interactions originated from the chargeable amino acid residues as well as the ability of the amino acid residues to form hydrogen bonds, causes peptide amphiphiles to self-assemble into a variety of nanostructures such as fibers, ribbons, bilayers, and micelles (Hendricks et al., 2017).

The self-assembly of peptide amphiphile (PA) molecules is controlled by the alkyl chain length, the formation of β-sheet hydrogen bonds and headgroup size and charge (Dehsorkhi et al., 2014). While the influence of the alkyl chain length and headgroup size on the self-assemble is relatively well understood for simple surfactants (Israelachvili, 2011), a fundamental understanding and control of the charge found in PA nanostructures are still lacking, particularly when the charge is not fixed but originates from amino acids with variable ionization states like glutamic acid (Deiss-Yehiely et al., 2017; Ortony et al., 2017). However, determining and modeling the effective charge density is crucial for understanding the PA self-assembly behavior and to optimize the sample processing conditions for a particular application. Furthermore, recent studies have shown that the molecular dynamics is critical for the PA nanostructures to form hierarchical structures (Freeman et al., 2018; Stupp et al., 2020; Wester et al., 2020) and to optimize biological signaling (Álvarez et al., 2021).

The filamentous peptide amphiphile nanostructures have been shown to mimic the extracellular matrix and have shown great potential for drug delivery and the regeneration of many different tissues (Matson and Stupp, 2012; Sato et al., 2018). Beyond their application as biomaterials, peptide amphiphiles have also been explored in material science applications and have been used to construct anisotropic actuating materials (Chin et al., 2018; Li et al., 2020). Consequently, there is an increasing number of experimental as well as theoretical studies into properties, behavior, and applications of peptide amphiphiles.

Experimentally, it is difficult to directly measure the amount of charge found in self-assembled peptide amphiphiles, since experimental techniques rely on indirect measurements. For example, ζ-potential measurements, a common tool employed in colloid chemistry, measures the electrophoretic mobility which can be used to estimate the electrostatic surface potential and charge. Also this method is best used for spherical particles and can it be more challenging to interpret for filamentous/fibrous structures. Ion counting or inductively coupled plasma mass spectroscopy (ICPMS) experiments (Gebala et al., 2015), can determine the number of condensed counterions and counterions contained in the cloud that surrounds a charged molecule or particle. ζ-Potential and ICMPS measurements provide an indirect indication of the strength of the electrostatic potential and charge of the charged moiety.

Small-angle X-ray scattering (SAXS) has been used to study supramolecular nanofibers (Bitton et al., 2005; Moyer et al., 2014; Kornmueller et al., 2016; Hendricks et al., 2017; Kornmueller et al., 2017). Anomalous small angle X-ray scattering (ASAXS) has emerged as a useful tool to determine the spatial distribution of ions surrounding charged macromolecules or self-assembled aggregates (Ballauff and Jusufi, 2006; Sztucki et al., 2012). This is done by careful subtraction of the scattered intensity profiles in the small-angle region at energies near and far away from an absorption edge of a targeted element such as rubidium or bromine. This technique has been used to determine the counterion distribution surrounding DNA (Das et al., 2003; Andresen et al., 2004), DNA-capped proteins (Krishnamoorthy et al., 2018) and spherical and cylindrical nanoparticles (Dingenouts et al., 2004; Patel et al., 2004).

Titration offers another method to establish (indirectly) the charge of self-assembling peptide amphiphiles. There are few titration studies of peptide amphiphile fibers and micelles (Toksoz et al., 2011; Buettner et al., 2017; Gao et al., 2017), due to bundling of the nanofibers at low or high pH values depending on whether the PA side chains are acidic or basic. Likewise, the possibility of morphological transitions make these titration experiments difficult. The interpretation of the titration curves is also complicated since a peptide amphiphile molecule can consists of multiple interacting amino acids. Furthermore the peptide amphiphile molecules are densely packed within the nanostructures and their ionizable side chains exist in a range of different chemical environments. Therefore (simple) relations, like the Henderson-Hasselbach (Atkins, 1998) equation, which estimates the amount of charge of a simple monoprotic acid in dilute solution, can not be applied to titration curves of self-assembled peptide amphiphiles.

Theoretically, molecular dynamics (MD) simulations have provided many insights into the self-assembly of peptide amphiphiles (Velichko et al., 2008; Lee et al., 2011, 2012). While MD simulations include many molecular details, they generally impose a fixed charge distribution and do not consider the possibility of dynamic chemical equilibria between protonation and deprotonation of acids. Acid-base chemical equilibrium can be introduced in MD simulations using constant-pH or reaction ensemble simulations (Donnini et al., 2011; Landsgesell et al., 2019), but this is difficult and costly to implement except in small, relatively simple systems (Donnini et al., 2011; Cote et al., 2014). Nanofibers, which involve many PA molecules, are therefore typically modeled with coarse-grained simulations, which impose fixed charges and other constraints and assumptions to predict the structure formation (Lee et al., 2011; Lee et al., 2012). Since self-assembling peptide amphiphiles often include amino acid side chains with weak acids like glutamic acid or aspartic acid (or weak bases like lysine), the extent of deprotonation or ionization is highly pH dependent. Another limitation is that MD simulations, due to their time-consuming nature, are not practical for systematic variation of parameters like salt concentration and pH.

Here, we present a theoretical model to describe the charge of a PA-nanofiber that includes the chemical equilibrium between the protonated, deprotonated, and ion-condensed states of the amino acid residue and the effect of dielectric environment and solvation. Importantly, the theory does not assume the charged state of the acid residues of the peptide but rather predicts the position-dependent state of charge. The theory is based on a molecular statistical thermodynamic approach that has previously been developed to predict thermodynamic and structural properties of end-tethered weak polyelectrolytes, which are polymers whose degree of deprotonation is not fixed but can change depending on environmental conditions including pH and ionic concentration (Solveyra et al., 2020). Predictions of this theory on the height of poly-acrylic acid brushes and the charged state of acid ligated gold nanoparticles were in good agreement with experimental observations (Gong et al., 2007; Wang et al., 2011). Likewise, a similar approach was used to investigate the effect of solution conditions on the charge regulation of bacteriophage capsids (Nap et al., 2014a) and the absorption of acidic ligands to quantum dots (Westmoreland et al., 2019). Also noteworthy is a recent study by Tagliazucchi et al. in which the molecular theory was used to study the self-assembly of neutral and chargeable peptide amphiphiles (Zaldivar et al., 2018; Zaldivar et al., 2019). However, ion-condensation and effects related to the electrostatic solvation free energy were not considered.

The objective of this theoretical study is to gain more insight into the effect that pH and other solution conditions, like type and amount of counterions, have on the amount of charge in the PA-nanofiber. To this end we consider a peptide amphiphiles with sequence C ₁₆ − V ₂ A ₂ E ₂, which is well characterized experimentally (Goldberger et al., 2011; Ortony et al., 2014).

One of the major findings is that there is a very strong coupling between the acid-base equilibrium and the local dielectric environment, which results in a much lower degree of deprotonation of the glutamic acid groups of the PA-nanofiber compared to the degree of deprotonation of the same acid in solution. Assuming a constant dielectric background, a common approximation employed in theoretical descriptions of colloids, polyelectrolyte solutions, and end-tethered polyelectrolyte layers, results in qualitatively and quantitatively very different degrees of deprotonation of the acid groups as compared to considering a position-dependent dielectric constant plus electrostatic solvation energy.

The paper is organized as follows. First, we present the theoretical approach. This is followed by a presentation and discussion of the results. Finally, we end with a summary, concluding remarks, and indicate potential future directions.

Here we consider C ₁₆ − V ₂ A ₂ E ₂, a well studied peptide amphiphile, consists of an aliphatic tail of 16 hydrocarbons linked to a peptide sequence of two valines, two alanines, and two glutamic acids, which forms a nanofiber. Figure 1 shows a “schematic” of the PA-nanofiber. The self-assembling PA-nanofiber is in contact with an aqueous solution. The reservoir is characterized by a given pH and contains either NaCl to mimic the high salt of biological media or RbCl since Rb⁺ provides a higher atomic number with better contrast for X-ray studies like anomalous small-angle X-ray scattering (ASAXS) (Ballauff and Jusufi, 2006). In this model, we assume complete dissociation of the salt ions. The pH is adjusted by adding either HCl or RbOH. In case of an aqueous solution that contains NaCl salt, either HCl or NaOH is added.

The carboxylic acid of the glutamic amino acid residue can be found either in a deprotonated state (A⁻), a protonated state (AH) or a state in which the acid is condensed with either a Rb⁺ or a Na⁺ counterion, which are denoted by ARb and ANa, respectively.

Following chemical reactions are included A H ⇌ A − + H + , (1) A − + N a + ⇌ A N a , (2) A − + R b + ⇌ A R b . (3) Considering both Rb⁺ counterions as well as Na⁺ ions allows us to investigate how different monovalent ions affect the protonation state of PAs.

The free energy describing one PA-nanofiber in contact with an electrolyte solution has a number of distinct contributions, which can be summarized as follows F = − T S conf,PA − T S mix,PA + F hb,PA + E VdW − T S mix + F chem + F elect + E elect, solv + E rep . (4)

The first four contributions are related to the conformational and translational entropy of the PA molecules, the van der Waals or hydrophobic interactions among the hydrocarbons and amino acids and the possibility to from hydrogen bonds between the amino-acids residues, respectively. The following four contributions are stemming from the mixing entropy of the mobile ions and solvent (S _mix), the acid-base chemical equilibrium of the acid group of the glutamic acid and ion condensation (F _chem), the electrostatic interaction energy (F _elect), and the electrostatic solvation energy of the charged species (F _{elect, solv}). The last term, E _rep accounts for the steric repulsions among all molecular species. We shall invoke a simplifying approximation. Therefore, we shall only present the pertinent terms and relegate the presentation of the complete theory to the supporting material.

The − TS _mix in the free energy corresponds to the mixing or translational entropy per unit length of the PA-nanofiber of the solvent (water) and mobile ionic species − S m i x k B L = ∑ k ∫ d r G r ρ k r ln ρ k r v w − 1 . (5) The index k runs over all the different types of mobile species: the water molecules, cations (Rb⁺, Na⁺), anions (Cl⁻), protons (H⁺), and hydroxide ions (OH⁻). The variable ρ _k (r) corresponds to the number density of mobile species k and v _w is the volume of a water molecule, which is used as the unit of volume. We use cylindrical coordinates to reflect the symmetry of the PA-nanofiber, and assume the system to be laterally homogeneous and only explicitly anisotropic in the radial direction r. The function G(r)dr is the cylindrical volume element divided by the unit length of the PA-nanofiber and equals G(r)dr = A(r)dr/L = 2πrdr. The next term in the free energy, F _chem, describes the chemical free energy associated with (de)protonation of the acid of the Glu amino acid and the ion condensation of Rb⁺ and Na⁺. F chem k B T L = ∫ d r G r ⟨ ρ G l u r ⟩ f A − r ln f A − r + β μ A − ⊖ + f A H r ln f A H r + β μ A H ⊖ + f A N a r ln f A N a r + β μ A N a ⊖ + f A R b r ln f A R b r + β μ A R b ⊖ + ∑ k ∈ H + , O H − , R b + , N a + β μ k ⊖ ∫ d r G r ρ k r . (6) Here, the f A − ( r ) is the fraction of glutamic acid residues that are deprotonated at position r, f _AH (r) is the fraction of protonated acids, and f _ARb (r) and f _ANa (r) are the fractions of acids that are condensed with either Rb⁺ or Na⁺. The variable μ i ⊖ corresponds to the standard chemical potential for a molecule of type i. The first and third terms within Eq. 6 describe the entropy of the deprotonated charged state (A ⁻) and protonated state (AH), respectively. The second and fourth terms in Eq. 6 correspond to the standard chemical potential of the charged and uncharged state, respectively. The subsequent terms have an identical meaning and pertain to the enthalpic and entropic contribution arising from the ion-condensation of Rb⁺ and Na⁺.

The seventh term, F _elect, in the free energy functional, Eq. 4 describes the electrostatic contribution to the free energy (Schwinger et al., 1998) and is given by F elect L = ∫ d r G r ⟨ ρ q r ⟩ ψ r + 1 2 ϵ 0 ϵ r r ∇ r ψ r 2 . (7)

Here, ψ(r) is the electrostatic potential and ⟨ρ _q (r)⟩ is the total charge density. The total charge density is the sum of the charge density of all charged mobile ions and the charge density of the amino-acid residues 〈 ρ q r 〉 = ∑ i e z i ρ i r + e z A − f A − r 〈 ρ G l u r 〉 . (8)

Here, the summation runs over all charged mobile ions Rb⁺, Na⁺, Cl⁻, H⁺, OH⁻ with z _i corresponding to their valency. e is the unit of charge. The second term is the number density of amino-acid residues that are deprotonated: ⟨ ρ A − ( r ) ⟩ = f A − ( r ) ⟨ ρ G l u ( r ) ⟩ .

In the electrostatic functional, ϵ ₀ and ϵ _r (r) correspond to the dielectric permittivity of vacuum and the local relative dielectric constant, respectively. Variation of the above functional with respect to the electrostatic potential yields the Poisson equation (Schwinger et al., 1998; Wang, 2008b). As shown by Z.-G.Wang, using field-theoretical arguments, (Wang, 2008a, 2010), consideration of a varying dielectric constant, that is having a non-uniform inhomogeneous dielectric medium, also requires the inclusion of a non-uniform electrostatic self-energy into the free energy functional. This electrostatic solvation energy (Wang and Wang, 2014) is the electrostatic energy required to transfer a charged molecule into a dielectric medium with a given position-dependent dielectric constant and is represented through a Born-type solvation energy. It is given by E elect, solv L = ∑ k ∫ d r G r Δ u k B r ρ k r with Δ u k B r = z k 2 e 2 8 π ϵ 0 a k 1 ϵ r r − 1 ϵ w . (9) For every ion of type k, a _k corresponds to the radius and z _k represents its valence. In order to describe the dielectric medium we need to provide a constitutive equation for the relative dielectric constant. We assume that the relative dielectric function of the medium is the volume-weighted average of the dielectric constant of water and the PA-nanofiber. ϵ r r = ∑ i ϵ i ρ i r ≈ ϵ P A 〈 ϕ P A r 〉 + ϵ w 1 − 〈 ϕ P A r 〉 (10) Here, we ignore the effect that the local volume fractions of the ions have. The dielectric constant of water equals ϵ _w = 78.54, while we set the dielectric constant of the PA-fiber equal to ϵ _PA = 2 (oil).

Using a uniform dielectric background media with a dielectric constant equal to that of the solvent is a common approximation in classical density functional theories (DFT) and Poisson-Boltzmann approaches that describe electrolyte solutions or electrolyte solutions near interphases, because the density of the electrolytes is low (Avni et al., 2019). A uniform dielectric background media is also employed in modeling end-tethered polyelectrolyte layers and is shown to work reasonably well, even for polyelectrolyte layers with higher grafting densities (Nap et al., 2006; Nap et al., 2014b). Electrostatic solvation free energy contribution, albeit constant, also occurs in, for example, the modeling of electrolyte adsorption to liquid-liquid interphases (Kung et al., 2009). Considering that the PA-nanofiber has a very dense aliphatic core assuming a uniform dielectric constant might be questioned. Therefore, we shall consider the case of a varying dielectric constant as well as a uniform dielectric background.

The last term in Eq. 4, E _rep, describe the repulsive interactions between all molecules, which, in the theory, are modeled as excluded volume interactions. These intermolecular excluded volume interactions are accounted for by assuming that the system is incompressible at every position ⟨ ϕ P A r ⟩ + ∑ k ϕ k r = 1 . (11) Here ϕ _k (r) = ρ _k (r)v _k is the volume fraction of mobile species k with volume v _k and ⟨ϕ _PA (r)⟩ is the volume fraction of the PA-fiber. These volume packing constraints are enforced through the introduction of the Lagrange multipliers π(r) since these are constraints, they are not formally part of the Helmholtz free energy. The volume fraction of the PA-nanofiber is the volume weighted average of the density of all residues of the peptide amphiphiles 〈 ϕ P A r 〉 = ∑ i ∈ A − , A H , A R b , A N a 〈 ρ G l u r 〉 f i r v i + ∑ r ρ r r v r . (12)

The first term describes the volume associated with the glutamic acid residues. Since we allow for volume changes upon ion-condensation to the glutamic acids, a summation over all charged states of the glutamic acid residues is required. The second contribution in Eq. 12 represents all other residues and the hydrocarbon tail.

The total free energy is minimized with respect to the number density of all species, ρ _i (r), the fraction of the charged states, f _k (r), and varied with respect to the electrostatic potential, ψ(r), under the constraints of incompressibility and the fact that the system is in contact with a bath of cations, anions, protons, and hydroxide ions. Thus, the proper thermodynamic potential is the semi-grand potential (Nap et al., 2006; Nap et al., 2014b): Ω = F − ∑ _i μ _i N _i , with N _i denoting the total number of particles of species i.

It is important to note that the complete theory, as formulated in the supporting material, is also dependent on the volume fraction of the PA-nanofiber, ⟨ϕ _PA (r)⟩, and the Glu number density, ⟨ρ _Glu (r)⟩. Meaning the total free energy includes the conformational energy and van der Waals interactions of the PA as well. However, due to technical challenges, we have not been unable to solve the complete molecular theory for appropriate large PA aggregation numbers yet. At these conditions the solvent is expelled from the aliphatic core and it becomes numerical challenging to satisfy the packing constraint. Hence we invoked an additional approximation. Namely, we imposed the volume fraction distribution of the PA-nanofiber and the distribution of the Glu number density. Similar approximations have been employed to model the counterion distribution around Au-NPs coated with DNA (Kewalramani et al., 2013) and as well as protein-nucleic acid conjugates (Krishnamoorthy et al., 2018).

To obtain a reasonable estimate of the volume fraction of the PA-nanofiber and the Glu number density, we use predictions from atomistic MD simulations. Technically, the simulations provide the number density of the water ( ρ w 0 ( r ) ) from which we can obtain the volume fraction of the PA-nanofiber needed as input for the theory. Ignoring the volume contributions from the ions, the volume fraction of the PA-nanofiber is then simply given as ⟨ ϕ P A 0 ( r ) ⟩ = 1 − ρ w 0 ( r ) v w . Because we allow for volume changes in the ion condensation reactions, we take the volume fraction of the PA-nanofiber from the simulation to correspond to the deprotonated state of glutamic acid. Consequently the total volume fraction of the PA-nanofiber is given by ⟨ ϕ P A r ⟩ = ⟨ ϕ P A 0 r ⟩ + ∑ i ∈ A H , A R b , A N a ⟨ ρ G l u r ⟩ f i r v i − v A − . (13)

Given the volume fraction profile and distribution of glutamic acid residues, we can compute using our theory the amount of charge of the Glu residues on the PA-nanofiber as well as the distribution of ions and water molecules and the electrostatic potential.

Finally, there has been considerable amount of theoretical research on charge regulation on flat surfaces (Ninham and Parsegian, 1971; Behrens and Borkovec, 1999; Boon and van Roij, 2011) and colloidal particles (Carnie et al., 1994; Stankovich and Carnie, 1996; Pericet-Camara et al., 2004; Frydel et al., 2007; Popa et al., 2010a; Popa et al., 2010b; Markovich et al., 2017; Avni et al., 2019). These are mostly formulated within the context of Poisson–Boltzmann theory, where the surface charge is either imposed to obey the acid-base equilibrium based on its solution equilibrium constant or by invoking the so-called the constant regulation approximation, (Behrens and Borkovec, 1999; Pericet-Camara et al., 2004; Popa et al., 2010a), that involves an unknown adjustable parameter. The approach presented here is, however, different. First, we consider a distribution of chargeable sites, i.e., there is a “diffusive” interface instead of a flat, impenetrable surface. Secondly, we do not impose the chemical state instead it follows from the free energy minimization. Thirdly, we considered a varying dielectric constant and electrostatic solvation energy. Common Poisson-Boltzmann based theories usually assume a fixed dielectric constant and only involve the translational entropy of the ions (F _mix) and electrostatic interactions (F _elect). The other terms for the free energy are usually not considered, especially the ones describing the chemical equilibrium of the acid groups, the ion-condensation of the ions and the effect of the dielectric environment as well as the excluded volume repulsion or finite size of all molecules-including water. Thus, unlike traditional Poisson-Boltzmann approaches, the current model includes the effect of charge regulation and dielectric environment as well as the type and size of ions.

Minimization of the free energy yields the following expression for the local volume fraction of the solvent ϕ w r = ρ w r v w = exp − β π r v w , (14) while the density of the ions reads ρ γ r = 1 v w exp β μ γ − β μ γ ⊖ − β π r v γ − β ψ r z γ e exp − β Δ u k B r . (15)

Observe that the Lagrange multiplier, π(r), can be interpreted as the lateral osmotic pressure. Also, notice that the chemical potential of water is not specified explicitly because the incompressibility constraint reduces the number of independent thermodynamic variables. Therefore, the chemical potential, μ _γ , is in reality an exchange chemical potential, i.e., the difference between the chemical potential of the molecule of type γ and that of water. Likewise, the charge neutrality and the water self-dissociation equilibrium further reduce the number of independent thermodynamic variables. Values of the exchange chemical potentials for all species can be obtained by relating them to their reservoir concentrations: ρ γ b u l k v w = exp ( β μ γ ⊖ − β μ γ − β π b u l k v γ ) , with Δ u k B , b u l k = 0 (Eq. 9) and ψ ^bulk = 0 (Eq. 16) (Nap et al., 2006; Nap et al., 2017).

Functional variation of the free energy with respect to the electrostatic potential yields the Poisson equation and its boundary conditions − ϵ 0 1 r d d r ϵ r r r d ψ r d r = ⟨ ρ q r ⟩ ; − ϵ 0 ϵ r r d ψ r d r r = 0 = 0 ; lim r → ∞ ψ r = 0 . (16)

Minimization of the free energy with respect to the different states of the Glu residue, f A − ( r ) , f _AH (r), f _ANa (r), and f _ARb (r) results in following set of ‘chemical reaction’ equations f A − r f A H r = e − β Δ G A H ⊖ e − β Δ G A H s o l v r e − β π r Δ v A H ρ H + r v w , (17) f A − r f A N a r = e − β Δ G A N a ⊖ e − β Δ G A N a s o l v r e − β π r Δ v A N a ρ N a + r v w , (18) f A − r f A R b r = e − β Δ G A R b ⊖ e − β Δ G A R b s o l v r e − β π r Δ v A R b ρ R b + r v w . (19)

The variable Δ G i ⊖ is the standard reaction free energy change of either the acid-base equilibrium reaction of the acid or the dissociation equilibrium reaction of the metal-ion pairs: ANa, or ARb. Here, Δv _i corresponds to the difference in volume between the products and reactants. Thus, Δ G ⊖ = ∑ J ν J μ J ⊖ and Δv _i = ∑ _J ν _J v _J . Here ν _J is the stoichiometric coefficient of species J involved in either the acid-base or ion-condensation reaction. The variable β Δ G A H s o l ( r ) is a position-dependent “solvation” free energy change and is given by Δ G A H s o l r = Δ u A − B r + Δ u H + B r − Δ u A H B r + E s o l v r v A − − v A H (20)

Δ G A N a s o l v ( r ) and Δ G A R b s o l v ( r ) are defined in a similar fashion. Note, that the expression for Δ G i s o l vanishes if a uniform dielectric background is assumed. In the above equation E _solv (r) is given as E s o l v r = 1 2 ϵ 0 ϵ r ′ ϕ P A r ∇ r ψ r 2 + ∑ k u k B r ρ k r ϵ r ′ ϕ P A r ϵ r ϕ P A r . (21) with ϵ r ′ [ ϕ P A ( r ) ] corresponding to the functional derivative with respect to the volume fraction of the PA-nanofiber and u k B ( r ) = z k 2 e 2 / 8 π ϵ 0 ϵ r ( r ) a k , i.e., to “absolute” electrostatic solvation energy. In our previously studies of weak polyelectrolyte layers, this solvation energy contribution did not appear since we only had considered the acid-base equilibrium, for which we assume v A − = v A H . The ion condensation has a contribution similar to Eq. 20, with a non-zero volume change equal to ( v A − − v A R b ) .

The Δ G i ⊖ is related to the chemical equilibrium constant, namely K i ⊖ = exp ( − β Δ G i ⊖ ) . Explicitly, the standard free energy change of the acid-base equilibrium AH ⇌ A ⁻ + H ⁺ is given by Δ G A H ⊖ = μ A − ⊖ + μ H + ⊖ − μ A H ⊖ , and the change in volume is equal to Δ v A H = v A − + v H + − v A H . The chemical equilibrium constant K A H ⊖ is related to the experimental acid-base equilibrium constant K A H = C ⁡ exp ( − β Δ G A H ⊖ ) of a single acidic monomer in infinitely dilute solution. Here C is a constant required for consistency of units and equal to C = 1/(N _A v _w ), where N _A is Avogadro’s number. This constant can be obtained readily by recasting the above chemical reaction into A − r H + r A H r = K A H ⊖ N A v w e − β Δ G A H s o l v r e − β π r Δ v A H , (22) and observing that, in the limit of infinite dilution π → 0 and ψ → 0, the two exponents vanish and the right-hand side becomes equal to the experimental equilibrium constant of a carboxylic acid monomer in solution. The chemical equilibrium constant for the ion condensation reactions are obtained in an identical fashion. Inspection of above reaction equations reveals that with increasing lateral osmotic pressure, π(r), corresponds to decreasing volume fraction of water or increasing peptide amphiphile volume fraction, and negative electrostatic potential the acid-base equilibrium shifts towards the neutral state. A similar trend occurs with decreasing dielectric constant (increasing ⟨ϕ _PA (r)⟩), in which the acid-base equilibrium shifts towards the neutral state.

To obtain a solution we need to solve the Poisson equation and incompressibility constraint simultaneously, since the unknowns in Eqs. 14, 15, 17, 18, 19 are the Lagrange multipliers or lateral pressures, π(r), and the electrostatic potential, ψ(r). Once the lateral pressures and the electrostatic potential are established, the amount of protonation and ion condensation of the Glu acids as well as the density distribution of the ions and the solvent are known. Solutions can be obtained by substituting the expressions of the volume fractions of all components into the incompressibility constraint Eq. 11 and the Poisson equation Eq. 16. This results in a set of non-linear integrodifferential equations that are converted into a set of coupled non-linear algebraic equations by discretizing space and these non-linear equations can be solved numerically (Hindmarsh et al., 2005). For detailed discretized expressions for individual terms and the computational procedure, the reader is referred to the supporting materials and previous publications (Nap et al., 2006; Nap et al., 2014a).

The inputs required to solve the non-linear equations are the concentrations of RbCl and NaCl, the reservoir pH the volume of all species, the volume fraction distribution of the PA-nanofiber and the distribution of chargeable sites of the Glu-residues. Also required are the acid-base equilibrium constants pK _a of the carboxylic acid and equilibrium dissociation constant of the Rb⁺ and Na⁺ ions with the carboxylate groups, all listed in the supporting information.

Figure 1 shows a representative snapshot of the atomistic MD simulation of the peptide amphiphile self-assembled into a nanofiber structure. Here the alkyl chains comprise the core of the nanofiber and the water-exposed periphery is made up of the peptides. For clarify we also show the chemical structure of the peptide amphiphile. Atomistic molecular dynamics simulations were performed based on the previously equilibrated structures (Ortony et al., 2017). The simulation details are provided in the Supporting Materials. Shown in Figure 1C are the volume fraction distribution of the PA-fiber and the number density of the Glu units.

Using the above volume fraction profile and distribution of glutamic acid residues as input, we determine with the theory the amount of charge on the PA-nanofiber, the distribution of ions and water molecules and the electrostatic potential. To characterize the charging behavior of the PA-nanofiber and quantify the changes that occur due to changes in pH and salt concentration, we computed the average degree of dissociation. The average degree of dissociation ⟨f⟩ is obtained by integration of the local position-dependent degree of dissociation, f A − ( r ) (Eq. 17), and given by ⟨ f ⟩ = ∫ d r G r f A − r ρ G l u r / ∫ d r G r ρ G l u r . (23)

Figures 2A,B show the average degree of charge as a function of pH for PA-nanofiber immersed in aqueous solutions that have different concentrations of RbCl and NaCl, respectively. For comparison, the dotted line labelled ‘ideal’ shows the degree of charge of an isolated carboxylic acid molecule in dilute solution, which obeys the ideal solution chemical equilibrium equation: ⟨ f ⟩ = 1 / ( 1 + 1 0 ( p K a − p H ) ) .

There is a very large deviation from ideal solution behavior as well as a strong salt dependence. For all salt concentrations and pH values, the degree of dissociation of the acidic residues is significantly lower than predicted by ideal solution behavior. For example at pH = 7 a carboxylic acid in solution is almost completely charged (⟨f⟩ = 0.99), since its pK _a equals 5. But the same acid residue in the PA-nanofiber has a degree of deprotonation reduced to less than 50%. To be precise, ⟨f⟩ = 0.56 for [RbCl] = 200 mM and by decreasing the salt concentration to 1 mM the fraction of charged acids further reduces and becomes as low as 15.3%. Thus, the degree of charge as a function of pH or the titration curve shifts to much higher pH values. The apparent pK _a , that is the pH for which half of the acid units are charged is shifted upwards by as much as 2–4 pH units, depending on the salt concentration. Also, we observe the occurrence of a plateau value of ⟨f⟩ = 0.8 at higher pH values. The fraction of charged acids never reaches 1. Thus, the amino acids inside a PA-nanofiber behave quantitatively and qualitatively differently from an isolated acid group.

To explain the predicted trends, we first focus upon the case of relatively low concentration of [RbCl] = 50 mM and present in Figure 3 the average fraction of acids that are protonated, deprotonated, and condensed with Rb⁺ as a function of pH. At low pH, the acids are mostly found in their protonated state. With increasing pH, the acid starts to gradually deprotonate. Observe that significant deprotonation only occurs for pH values well above the pK _a of the carboxylic acid group. Simultaneously with deprotonation ions start to condense. The general behavior of the degree of dissociation as a function of salt concentration and pH can be understood as follows. The charge or degree of deprotonation of the amino acids occurs through a balance between the chemical free energy of the acid-base equilibrium reaction and the counterion ion condensation to the glutamic amino acids, the electrostatic interactions, and the mixing entropy of the solvent and mobile ions. To be exact, the entropy that is associated with counterion confinement. With decreasing salt concentrations, the electrostatic repulsion between the deprotonated amino acids are less screened and the system experiences increased electrostatic repulsion, for which the system tries to compensate. The PA-nanofiber system can respond in several ways.

First, the PA system recruits additional counterions from the reservoir, which increases the electrostatic screening. This decreases the (enthalpic) electrostatic repulsion at an entropic cost related to the loss of translational entropy of the counterions. Secondly, the electrostatic repulsion can be reduced by decreasing the number of deprotonated amino acids, which can be accomplished by either shifting the acid-base equilibria towards its protonated state or by condensing counterions. Both occur at the cost of their respective chemical free energies. For a weak acid the cost associated with the chemical free energy contribution is small as compared to the loss of entropy due to counterion confinement. Hence the mechanism of charge regulation through shifting the chemical equilibrium is usually the primary mode by which charged systems, such as weak polyelectrolytes and PA-nanofibers, try to mitigate the effects of electrostatic repulsion.

There is a third mechanism available to the peptide-amphiphiles to negate the electrostatic repulsion by increasing the intermolecular spacing between PA molecules. By changing the spatial distribution of the PA the charges are moved apart and the electrostatic repulsion decreases. Clearly, the approximation we use does not account for this possibility, since we assume a fixed spatial distribution of Glu residues within the PA-nanofibers. Observe that increasing the intermolecular spacing between PA molecules cannot negate nearest neighbor electrostatic interactions of charged Glu residues that reside on the same PA molecule. Also, the hydrophobic interaction of alkyl chains leads to the formation of the very dense core of the PA-nanofiber, this limits the possibility of increasing the intermolecular spacing between PA molecules within a (cylindrical) nanofiber. The possibility of molecular reorganization depends also on the strength of the β-sheet hydrogen bonds of the peptides. Only by a large structural transition, which involves enthalpic and entropic penalties, from a cylindrical nanofiber to a spherical micelle can the electrostatic penalty be reduced (Israelachvili, 2011; Zaldivar et al., 2018; Zaldivar et al., 2019). This occurs experimentally for pH above 8. We have investigated the effect of different spatial distributions of PA molecules on the charge of the PA-nanofiber. The distributions are obtained from MD simulations with different aggregation numbers. Calculations using different distributions resulted in qualitatively similar titration curves. We conclude that changes in spatial reorganization of the PA-nanofiber have a minor effect in determining the overall charge of the PA-nanofiber.

We have developed a theoretical model to describe and predict the charges found in self-assembled peptide amphiphile nanofibers as a function of pH and ion concentration. In particular, we computed the amount of charge on a nanofiber made of peptide amphiphiles with the sequence C ₁₆ − V ₂ A ₂ E ₂. Theoretically, we accounted for the acid-base equilibrium as well as the local dielectric environment by allowing for a position-dependent dielectric constant as well as a local solvation energy of the charged species. We find that the charge on the glutamic acid residues is much lower compared to the same acid in dilute solution. Surprisingly there is a very strong coupling between the acid-base equilibrium and the local dielectric environment. Considering a constant dielectric background instead of varying dielectric media results in qualitative very different degrees of deprotonation. Also, the shape and value of the electrostatic potential and counterion ion distribution are quantitatively and qualitatively different. Considering a constant dielectric media, a commonly applied approximation, result in counterion penetration into the alkyl region of the peptide amphiphile nanofiber. This erroneous result indicates that it is necessary to include a dielectric constant plus electrostatic solvation to properly describe the charging behavior of peptide amphiphiles.

We also considered the effect of counterion ion binding of Rb⁺ and Na⁺ to the carboxylate group of glutamic acid. We demonstrated that ion condensation can result in a considerable amount of bound ions, particularly for buried carboxylic groups. Counterion binding offers an additional mode of charge regulation and, depending on the strength of the metal-acid binding, a “bimodal” distribution of counterions is found with spatially distinct regions of bound and free of counterions. Experimentally, this counterion distribution could be revealed by X-ray scattering experiments.

It can potentially also be probed through nuclear magnetic resonance (NMR) measurements or inductively coupled plasma mass spectroscopy (ICPMS). NMR has been used to measure the macromolecule-bound metal ion concentrations in for example protein and DNA metal complexes (Kozlyuk et al., 2016) while ICPMS has been used to reveal the concentration of counterions located around double-stranded DNA (Gebala et al., 2015). Observe that the latter method measure the total amount of counterions and can not distinguish between condensed and free counterions.

While the theory presented here includes many important details related to charge regulation and electrostatic solvation is still a mean-field approach and as such does not include electrostatic fluctuations. Another approximation is that we imposed the PA volume fraction and glutamic acid density distribution. Similar approximations have been used to investigate DNA-coated nanoparticles (Kewalramani et al., 2013). This assumption can be rationalized by noting that the PA-nanofiber structure is compact and dense, which is particularly true for the aliphatic core of the fiber. Because of this approximation, we can only predict the charge of cylindrical PA-nanofibers and cannot investigate the transition to for example spherical micelles. Experimentally such transitions are known to occur for elevated pH values, roughly above pH 8. Hence, the predictions are only representative for cylindrical nanofibers up to pH 8. We also investigated the charge of PA-nanofibers that had a different PA-volume fraction and Glutamic number density distribution. The distributions are obtained from MD simulations of PA-nanofibers with different aggregation numbers, which provides qualitatively similar results for the charging behavior of the PA-nanofiber. Therefore, we believe that the trends presented here are relevant and that changes in the PA distribution will have only a minor effect on the charging of PA-nanofibers. However, it is important to emphasize that because the model assumes a static peptide nanofiber, the results presented here are qualitative. The PA distribution and the amount of charge on the glumatic acid residues are coupled together and can therefore influence each other. This interplay between structure and charge of the PA molecules is missing since the PA molecular distribution is assumed to be static and cannot respond to changes in the charge distribution. For future directions, we will add the conformational entropy of the PA-molecules into the theory and investigate how the charging of PA-molecules couples with the self-assemble of the PA-molecules.

Another limitation of the theoretical approach is that the local dielectric constant is taken to be the linear volume-weighted average of the dielectric constant of water and the PA-nanofiber (Eq. 10). Such a constitutive equation provides an empirical description of the dielectric properties of the system. However, the local dielectric constant is coupled with the local polarization of the molecules, which is in turn determined by the local density composition and the local electric field. The description of the dielectric properties of the PA nanofiber could be improved upon by considering a polarizable molecular model. An example of such a model can be found in reference (Nakamura et al., 2012).

Very noteworthy are the recent theoretical investigations by Zaldivar et al., who applied the same Molecular Theory approach to self-assemble of peptide amphiphiles and demonstrated the structural transitions between bilayer, cylindrical, and spherical structured nanofibers can occur as a function of pH (Zaldivar et al., 2018; Zaldivar et al., 2019). Observe they opted to represent the peptide amphiphiles with much more coarse-graining. This much more coarse-grained representation, based up the Martini-force field, of the peptide amphiphiles enabled the successful solving of the equations. Another important difference was that their approach did not include the possibility of either a varying dielectric constant or position-dependent solvation energy. However, our results demonstrate that the inclusion of varying dielectric constant and solvation energy leads to significant changes in the charge of cylindrical peptide amphiphiles nanofibers, which should influence the location or “phase” boundary of structural transitions. How the changes in charge of PA-nanofiber influence the location and the possibility of structural changes is difficult to access, since other nanostructures like spherical micelles will also undergo a similar charge regulation mechanism. The calculations presented here demonstrate that the amount of charge and thereby the electrostatic potential is sensitive to the dielectric environment and this effect needs to be considered.

Finally, note that the all-atom simulations we used to obtain an input of PA-nanofiber distribution assume that both glutamic acid residues are completely charged and are independent of environmental conditions like pH and salt concentration. The results from our theory here show that the charge of PA-nanofibers is not fixed and varies considerably, depending on the surrounding environment. This suggest that all-atom simulations should utilize grand-canonical reaction ensemble (Landsgesell et al., 2019) or constant pH ensemble (Donnini et al., 2011; Cote et al., 2014) for an more accurate description of the charging of the PA molecules, instead of conventional all-atom simulations that assume a fixed charge.

The results, here, indicate that to properly understand the charge on peptide amphiphiles we need to take into account both the capacity of charge regulation of the amino acid residues and its coupling with the local dielectric environment.