Cross-linker effect on solute adsorption in swollen thermoresponsive polymer networks

The selective solute partitioning within a polymeric network is of key importance to applications in which controlled release or uptake of solutes in a responsive hydrogel is required. In this work we investigate the impact of cross-links on solute adsorption in a swollen polymer network by means of all-atom, explicit-water molecular dynamics simulations. We focus on a representative network subunit consisting of poly($N$-isopropylacrylamide) (PNIPAM) and $N$,$N'$-methylenebisacrylamide (BIS/MBA) cross-linker types. Our studied system consists of one BIS-linker with four atactic PNIPAM chains attached in a tetrahedral geometry. The adsorption of several representative solutes of different polarity in the low concentration limit at the linker region is examined. We subdivide the solute adsorption regions and distinguish between contributions stemming from polymer chains and cross-link parts. In comparison to a single polymer chain, we observe that the adsorption of the solutes to the cross-link region can significantly differ, with details depending on the specific compounds' size and polarity. In particular, for solutes that have already a relatively large affinity to PNIPAM chains the dense cross-link region (where many-body attractions are at play) amplifies the local adsorption by an order of magnitude. We also find that the cross-link region can serve as a seed for the aggregation of mutually attractive solutes at higher solute concentrations. Utilizing the microscopic adsorption coefficients in a mean-field model of an idealized macroscopic polymer network, we extrapolate these results to the global solute partitioning in a swollen hydrogel and predict that these adsorption features may lead to non-monotonic partition ratios as a function of the cross-link density.


Introduction
Responsive polymers have increasingly gained attention in many research fields due to their ability to reversibly adapt to external stimuli such as temperature, osmotic pressure, or pH.Various types and shapes at different length scales have been designed, providing various possibilities for applications, 1,2 such as solute uptake, transport 3 and release, 4,5 sensors, 6 intelligent coatings, switchable catalysis, 7,8 etc.To structurally stabilize and fine-tune properties and function, polymer architectures are often E-mail: joachim.dzubiella@physik.uni-freiburg.deequipped with chemical cross-linkers covalently interconnecting the chains, 9 which are then referred to as hydrogels.Typical responsive hydrogels in most studies barely exceed a molar cross-linker density of twenty percent, since greater values increase the rigidity of the gel and reduce the swelling properties due to the denser polymer network structure. 10,8 the zoo of constituents, thermoresponsive hydrogels based on poly(N -isopropylacrylamide) (PNIPAM) are among the most intensively investigated systems, since their volume phase transition at about room temperature as well as a high water content promise good biocompatibility 11,5 and make them convenient to handle.Pure PNIPAM was found to have the lower critical solution temperature (LCST) at roughly 304 K as reported by Heskins and Guilett in 1968. 12A frequently utilized cross-linker for 1 arXiv:1903.12043v2[cond-mat.soft]3 Apr 2019 PNIPAM gels, used in radical polymerization, is N ,Nmethylenebisacrylamide, often abbreviated as BIS or MBA.BIS has chemical similarity to PNIPAM (compare Fig. 1), is non-degradable, has a very high reactivity, and retains PNIPAM's LCST. 13,14,15,16,17esides these morphological properties, the degree of cross-linking influences solute uptake and partitioning.The partition ratio is the ratio of the solute concentrations inside and outside the gel and is therefore a crucial parameter controlling device functionality especially for drug delivery or catalytic systems.For the latter, for instance, metal nanoparticles inside hydrogels catalyze reactions and the effective reaction rates depend crucially on the concentration of the reactants in the permeable polymer matrix. 18,8,7,19he partition ratio may be affected by generic as well as specific cross-linker effects.The cross-linker density first of all simply changes the packing fraction and with that the overall steric exclusion by the polymer mesh. 11Furthermore, it has become clear that more complex, e.g., local attractive and/or electrostatic interactions can lead to complex and even cooperative effects in the partitioning. 20,21,22,23n particular, a 'vertex trapping' effect due to many-body attractions in the dense cross-link region has been reported in generic coarse-grained simulations of polymer networks. 24,25,26,27,28More specific chemical effects should also play a role, as indicated by all-atom molecular dynamics (MD) computer simulations of bare PNIPAM chains 29,30,31 , peptide-like chains 32,33,34 , these in combination with solutes with various polarity 35,36,37,38,39,40 as well as by simulations revealing the influence of cross-links to polymer networks solvation and structural properties. 41,42,43,44he aim of this work is to investigate the effects of crosslinking on solute adsorption in swollen hydrogels made up of PNIPAM and BIS (below the volume phase transition temperature (VPTT)) by utilizing all-atom, explicit water MD simulations of a minimal polymer network setup.In order to do this, we consider one BIS cross-linker with four atactic PNIPAM chains restrained in a tetrahedral geometry.In our analysis, we subdivide the solute adsorption regions and systematically distinguish between contributions stemming from polymer chains and cross-linker parts.We probe solutes of various polarity, representing typical chemical compounds, in the highly diluted regime.We finally demonstrate in a simple model, how these contributions affect the global solute partitioning in large hydrogels as accessible by experiments.

Hydrogel building blocks: PNIPAM and BIS
Constructing covalently cross-linked polymeric networks for our computer simulations requires two types of building blocks: chain monomers and the cross-linker (see Fig. 1).
The former provides two binding sites and builds up the chains.
The latter has four binding The partition coefficient may be affected by generic as well as specific crosslinker effects.The crosslinker density first of all simply changes the packing fraction and with that the overall steric exclusion by the polymer mesh. 11Further, it has become clear that more complex, e.g.local attractive and/or electrostatic interactions can lead to complex and even cooperative effects in the partitioning. 20,21][35] The aim of this work is to investigate the effects of crosslinking on solute adsorption in swollen hydrogels (below the LCST) by utilizing all-atom, explicit water MD simulations of a minimal hydrogel setup.For this we consider one BIS-linker with four atactic PNIPAM chains attached in a tetrahedral geometry.In our analysis, we subdivide the solute adsorption regions and systematically distinguish between contributions stemming from polymer chain and crosslink parts.We probe solutes of various polarity, representing typical chemical compounds, in the highly diluted regime.We finally demonstrate in a simple model of a hydrogel, how these contributions affect the global solute partitioning in a large hydrogel as accessible by experiments.sites and thus interconnects four chains.We have chosen poly(N -isopropylacrylamide) (PNIPAM) and N ,Nmethylenebisacrylamide (BIS) for the chains and the crosslinker unit respectively (see Fig. 1).

Hydrogel building blocks: PNIPAM and BIS
We employ the OPLS-QM2 force field recently developed by Palivec et al. 31 for the PNIPAM monomers.Compared to the standard OPLS-AA 46 parameters, this force field features a reparametrization of the partial charges retrieved from ab-initio calculations and further manual fine-tuning to reproduce the experimental LCST of PNIPAM.Due to the chemical similarity of BIS and PNIPAM, we adopt the very same partial charges for the cross-linker.These were confirmed by our own quantum mechanical calculations using the Gaussian 09 software. 47More details on the force field parameters are provided in Appendix A.

Setup
The actual polymerization procedure in experiments is subjected to randomness leading to different possible network topologies.Looking at the very local structure inside a hydrogel, the following sources for inhomogeneities are possible: dangling chains, entangled chains, and loops. 48This work, however, focuses on generic subunits of a defect-free network architecture of a swollen hydrogel, 22,49 namely the cross-linker and its four adjacent chains.Our setup consists of one BIS-linker connected with four PNIPAM chains (Fig. 2a), each composed of 12 monomers, i.e., 48 monomers in total.The terminal backbone carbon atoms of each chain are position-restrained in order to sustain a tetrahedral structure.It further facilitates a clear analysis and the results can be to some extent generalized, which will be dis-  --), which has an edge length of 4.69 nm and a vertex-to-centroid distance of L = 2.87 nm.For the sake of clarity, water molecules are not shown.(c) Probe solutes.Atoms are color-coded by element, i.e., black carbon, red oxygen, blue nitrogen, and white hydrogen atoms.See Table 1 for the full simulation specifications.
cussed later in this work.
The corners of the tetrahedron, generated by the position restraints, are L = 2.87 nm distant to its centroid and the edge length accounts to 4.69 nm (Fig. 2b).We have chosen this size to ensure a relative chain stretch λ between 0.75 and 0.85, which is expected in swollen hydrogels. 38The relative chain stretch λ is defined as the ratio of the mean end-to-end distance and the contour length The mean end-to-end distance of the chains is ee = 2.65 nm, where the brackets .. denote the ensemble average.The distance ee was measured from the cross-linker contact backbone atom to the chain's terminal backbone atoms during the N pT -simulation.The contour length per monomer ∆L c in an atactic PNIPAM chain is approximately 0.265 nm, and when multiplied by 12, one obtains the contour length of a single chain.Eventually, the average relative stretch in all simulations is λ = 0.83.Generating such a setup starts by placing the BIS-linker in the center of the box, which is the centroid of the virtual tetrahedron.At each backbone binding site, indicated by dashed lines in Fig. 1, PNIPAM monomers with random tacticity are attached in a head-to-tail manner with the backbone axis pointing towards the desired coordinates of the position restraint.PNIPAM's backbone bonds are squeezed to match the size of the tetrahedron and were allowed to relax during the first steps (energy minimization and equilibration) of the simulation.
The edge lengths of the rectangular simulation box are chosen large enough (6.32 nm × 7.65 nm × 9.64 nm on average) in order to ensure that the position-restrained back-bone terminals of different periodic images are separated by at least 3 nm in x-and y-, and 5 nm in z-direction (Fig. 2b).Thus, we avoid interactions between chain ends across the box boundaries and can locate a (water/solute) bulk phase in peripheral box regions.
The solutes are subsequently inserted at random positions in the simulation box, which is finally filled with more than 15000 SPC/E 50 water molecules.Details are listed in Table 1.

Probe molecules
We analyze the adsorption properties of several types of molecules covering different sizes and polarity.We focus on some aromatic compounds due to their application as model reactants in catalytic experiments 18,51 and since aromatic rings are found in numerous drugs, 52,53,54 e.g., common painkillers.Precisely, these compounds are benzene (B), nitrobenzene (NB), (uncharged) 4-nitrophenol (NP 0 ), and (charged) deprotonated 4-nitrophenolate (NP − ).
We further probe two alkanes, namely hexane (C 6 ) and butane (C 4 ), sodium chloride (Na + /Cl − ), and methanol (C 1 OH).All compounds are visualized in Fig. 2c.If not stated otherwise, we insert one probe molecule into the system to analyze the infinite dilution limit.To estimate finite concentration effects of aromatic compounds we perform simulations with 20 solutes.List of solutes and the simulation setups are listed in the summary Table 1.The standard OPLS-AA 46 force field was utilized except for the charged nitrophenolate NP − , for which the excess charge was distributed among the molecule due to the mesomeric effect, 38 leading to higher polarity of the nitro group.

Simulation details
We employed all-atom, explicit-water molecular dynamics (MD) simulations to study the polymer-solute interactions.The simulations were performed using the Gromacs 5.1 software 55,56,57 utilizing force fields mentioned above.
All covalent bonds of hydrogens were constrained with the LINCS 58 algorithm.The cut-off distance for Lennard-Jones and short-range electrostatic interactions was set to 1.0 nm while long range electrostatics was accounted for by the Particle Mesh Ewald (PME) method with cubic interpolation and a grid spacing of 0.12 nm. 59eriodic boundary conditions in all three directions were used and the simulations were carried out under constant temperature and pressure, which were controlled by the velocity-rescale thermostat (at T = 290 K, τ T = 0.1 ps) and the Berendsen barostat (at p = 1 bar, τ p = 1 ps), respectively. 60,61fter the initial energy minimization (steepest descent), the system was equilibrated in the N V T ensemble for 2 ns and in the N pT ensemble for another 10 ns.The integration step of the leap-frog integrator was set to 2 fs and data were collected every 10 ps.The total simulation time t sim per solute is summarized in Table 1.

Analysis and discussion
From the trajectories we calculate the center of mass (COM) positions of the cross-linker and the structure of the chain monomers, the solute(s) and water molecules around it.The distances r between the COM of the cross-linker and the COM of molecules is used to obtain the normalized radial density distributions g(r) = ρ(r)/ρ 0 for BIS-water and BISsolute.The bulk phase concentration ρ 0 (see Fig. 3 and Table 1) is obtained from simulations by calculating the number density along the z-axis and averaged in the region 1.5 nm distant to the restrained polymer atoms.
For the 20 NP − /Na + pairs, we only analyzed the nitrophenolate trajectories.In the case of the Na + /Cl − simulation, each ion type was analyzed individually.The results for sodium and chloride ions are very similar and yield the same results within the range of our precision and thus are presented for either type in Table 1.
Further, the (radial) PNIPAM monomer number distribution ρ mer (r) is retrieved, which helps us to distinguish between different polymer adsorption domains (Fig. 3).The solutes' distributions are the basis for calculating the solutepolymer adsorption in our setup as detailed below.We demonstrate how the splitting of the adsorption into the chain and the cross-linker contributions is achieved and how this can be used to estimate partition ratios of an entire hydrogel.

Polymer distribution
The position restraints restrict the movement of the PNI-PAM monomers and the cross-linker.The COM of the cross-linker in this setup fluctuates around the simulation box center with a mean displacement of ∆r xlink =0.23 nm.The PNIPAM distribution for each monomer has been evaluated with respect to the COM of BIS.The closest monomers to the cross-linker distribute in a bimodal fashion, stemming from multiple possible side-chain-side-chain and cross-linker-side-chain interactions.This effect averages out for chain monomers further distant from the linker resulting in Gaussian distributions.The distances of two adjacent monomers is roughly 0.25 nm.The distribution of all monomers together (i.e., averaged over a spherical shell ∝ ρ mer (r)4πr 2 dr including all chains) shows a plateau region between 1.5 and 2.5 nm, see Fig. 3b.In this range we find an almost constant monomer density, which can be used to evaluate the intrinsic adsorption per chain monomer.We thus define the number of monomers in an interval [r 1 , r 2 ] as with ρ mer (r) the radial density of PNIPAM monomers with reference to the cross-linker.The number has the upper bound N mer (0, ∞) = 48, which is the total number of PNI-PAM monomers in the system.

Solute adsorption in high dilution
The adsorption of solute inside the hydrogel of volume V gel (including the containing water) can be calculated from the The range of influence of the cross-linker is assumed to vanish around ra = 1.5 nm.Region II, ranging from ra to r b , is dominated by the linear PNIPAM chain.Here, the monomer number, ρmer(r)4πr 2 , is roughly constant.In interval III, from r b = 2.5 nm to rc = 3.5 nm, the polymer chains terminate.The average distance between the position-restrained backbone terminals and the cross-linker, L = 2.87 nm, is shown by the vertical dotted (red) line.For large distances, i.e., r > rc, we assume negligible influence from the polymer and consider region IV as 'bulk'.
Table 1 Simulation specifications and resulting solute adsorption coefficients to the polymer setup with a relative chain extension λ ≈ 0.83 as depicted in Fig. 2. Each N pT -simulation was carried out at T = 290 K, p = 1 bar, with more than 15000 water molecules, one solute (30 molecules and three pairs in the case of C 1 OH and Na + /Cl − , respectively) and analyzed within a simulation time of tsim.Nitro-aromatics were further tested with Nsolute = 20 molecules and we added Na + counterions for charged NP − .Note the differences in the simulation time, e.g., the simulations with Nsolute = 1 were carried out approximately ten times as long as with 20 individuals in order to reach sufficient sampling quality.Adsorption coefficients Γ * = Γ/ρ 0 are split into contributions as described by eqn ( 5) and ( 7) and are normalized by solute bulk concentration ρ 0 (Details are provided in the main text).For the NP − /Na + pairs, only results for NP − are presented.The ions Na + and Cl − were simulated together but analyzed separately.Both types yield very similar results.The most relevant results, the adsorption coefficients Γ * xlink and Γ * mer to the cross-linkers and monomers, respectively, are visualized in Fig. 5 Simulation Specifications Results radial distribution as

Aromatics
which is a Kirkwood-Buff integral 62,63 counting the excess (or deficit) number of solutes with respect to bulk concentration ρ 0 and depends on the volume V ex excluded by the polymer.The adsorption Γ = 0 refers to the scenario at which the attractive solute-polymer interaction fully compensates for the steric exclusion −ρ 0 V ex .
Transferring this concept to our setup (Fig. 2), for which radial density profiles g(r) of the solutes (Fig. 4a, b) are measured from the COM of BIS, we define the partial adsorption Γ(r 1 , r 2 ) counting excess solutes in the interval [r 1 , r 2 ], reading and can scan the adsorption in different domains with respect to the cross-linker as shown in Fig. 3.The total adsorption, i.e., Γ tot = Γ(0, ∞), is not only solute-specific but also depends on the number of monomers and the geometry.
To separate the effects of our particular system setup, we distinguish now between three different contributions, stemming from the cross-linker (Γ xlink ), linear chains (Γ chain ), and chain terminal ends (Γ end ).We will determine them by classifying different adsorption domains I, II, and III, and bulk phase (IV), as depicted in Fig. 3.The total adsorption can be written as In our setup, the total adsorption is dominated by the chain contributions due to the numerous PNIPAM monomers compared to only one cross-linker.The contribution of the chain ends Γ end is of lesser importance for this work.Dangling ends in hydrogels are very common but usually not of high concentration.It can be computed once the chain and cross-linker terms have been determined.In our setup, however, the calculated values of Γ end cannot be interpreted in a meaningful way due to the position restraints, which locally alter the relative water-polymer dynamics, in other words, disable the 'dangling' behavior of such terminals.
The equilibrium bulk concentration ρ 0 depends on the simulation box size and the binding affinity.It is convenient to define an infinite-dilution solute-specific adsorption coefficient that does not depend on concentration via The adsorption coefficients Γ * have the units of volume (nm 3 ), and correspond to highly diluted cases (ρ 0 → 0), in which solute-solute interactions can be neglected.In the case of the highly water-soluble methanol, tested with 30 molecules, and simple ions (three pairs), solute-solute interactions play a minor role for the adsorption.Thus they are to some extent considered as very diluted scenarios and are comparable to single-solute simulation results.The different adsorption coefficients for all compounds are summarized in Table 1.

Adsorption per chain monomer
The chain contribution to the solute adsorption is what one would expect from a single isolated linear PNIPAM chain, i.e., in the absence of the cross-linker and any other chains close by.It can be described by the adsorption per monomer Γ mer and with N mer (0, ∞) = 48 in our setup, this yields The adsorption per monomer is evaluated from the chain domain (where the BIS and end effects are negligible, see II in Fig. 3), i.e., r ∈ [r a , r b ], reading We now compare the adsorption of the solutes to the PNI-PAM chain, listed in Table 1 and visualized in Fig. 5.The results can further be compared with the density profiles (Fig. 4a, b).We start with the examination of the rather weakly adsorbing species (Fig. 4a).Methanol is the smallest probe molecule tested and is highly soluble in water and shows a rather low binding propensity.It is in fact slightly attracted to the polymer, but this cannot compensate the volume exclusion effect of PNIPAM and thus its adsorption coefficient is of negative value.Methanol's preferential adsorption has already been reported in experiments 64 and other simulations, 65,66,67,68 studying primarily the cononsolvency of PNIPAM in water-methanol mixtures.
Sodium and chloride have the lowest binding affinity to the hydrogel, which has already been shown in previous simulations of isolated chains. 38,35As expected, simple wellhydrated ions are repelled from low dielectric (less polar) regions.
The two probed alkanes, butane and hexane, have very similar profiles.The bigger hexane shows slightly higher adsorption than butane owing to the larger surface area, which facilitates hydrophobic interactions with apolar groups of the polymer chains.The very same argument does not hold when comparing with benzene.Benzene has about the same size as the alkanes, but shows higher binding affinity than the larger hexane.On the molecular level, the adsorption mechanism looks similar.Benzene and hexane tend to preferentially make contact with hydrophobic parts of the polymer.
Comparing the aromatic compounds, which are roughly of equal size, we find the adsorption generally increases with the polarity of their substituents.The order according to their polarity, starting with the apolar benzene, is B→NB→NP 0 →NP − , cf.Fig. 4b.All of them own a hydrophobic aromatic ring, interacting with the hydrogel described as in the benzene case.With one polar substituent for benzene, namely the nitrobenzene, the adsorption is more than doubled.This stems from additional hydrogen bonding 45 between the nitro-oxygens and the polymer's amide hydrogens.The very same interaction mechanism applies to NP 0 and NP − .The extra hydroxy tail in the case of nitrophenol (NP 0 ) does not lead to a significant change of the adsorption to PNIPAM.On the one hand, the OH group can interact with the polymer's amide group and on the other hand, increases the water solubility.These two effects seem to compensate for NP 0 adsorption to the chains, such that the adsorption is similar to the one for NB.
The deprotonated and hence charged NP − is the best adsorbing compound tested.The deprotonation leads to a redistribution of the electronic density, increasing the polarity of the whole molecule.The higher charging of the nitro-oxygens as well as the O − tail stabilize the contacts with PNIPAM's amide hydrogens, resulting in a roughly two times higher adsorption coefficient compared to NB and NP 0 .

Effects of the cross-linker
The contribution of the cross-linker in eqn ( 5) is obtained by integrating the solute radial density distribution from the COM of the BIS-linker up to the onset of the chain domain r a (cf.domain I in Fig. 3) and subtracting the estimated linear chain contribution therein, formally written as xlink in the infinite dilution regime, as summarized in Table 1.The adsorption to the monomers increases with size, e.g., compare C 4 → C 6 .The aromatic molecules show the highest adsorption, which is magnified by increasing polarity B → NB → NP 0 → NP − , coinciding with a crosslinker-enhanced binding affinity.
This quantity combines specific interactions of the solute with the BIS-linker and the more complex many-body effects resulting from the higher concentration of PNIPAM monomers.
Note that in comparison to PNIPAM monomers, BIS has two amide groups and no isopropyl groups, thus creating a more hydrophilic environment than the chains.The apolar compounds, B, C 4 , and C 6 show a slightly negative binding affinity in the cross-link region.In contrast, see again Table 1 and Fig. 5, the solute adsorption Γ * xlink increases with polarity, where nitro-aromatic solutes are especially attracted.Nitrobenzene shows a more than doubled adsorption to the cross-link region when compared to bare chain monomers.The NP 0 and NP − adsorption per cross-linker is even tenfold higher.The binding mechanism is similar to the single chain adsorption.The numerous amide hydrogen combinations make it very probable for the nitro-oxygens to find binding partners.
As already stated for the chain adsorption, NP − has the most polar nitro group resulting in the strongest adsorption coefficients in this study.Examining the simulation trajectories, we repeatedly found NP − in the location shown in Fig. 6.One or both of the nitro-oxygens couple (forming hydrogen bonds) with two to three hydrogens from the amide groups: one from BIS and one or two from the PNIPAM monomers.Additionally, the hydrophobic isopropyl groups or the backbone of PNIPAM can contact, almost embed, the aromatic ring, enhancing the stability of such an adsorbed state.The same mechanism has been observed for NB and NP 0 , but the higher partial charges of NP − promote the binding.

Finite concentration effects
The adsorption in the low density case can differ from scenarios with higher concentrations owing to solute-solute interactions.This was tested with nitrobenzene, nitrophenol, and nitrophenolate, using twenty molecules per species, where we moved by an order of magnitude from the 2-3 mM concentration range up to 30-40 mM.The concentrations and local adsorption results are also summarized in Table 1, while density profiles and structures are shown in Fig. 4b,  c, where we compare them to the low-density limit.We find for all tested solutes significant collective effects.The linear dependence of the adsorption on ρ 0 thus only holds for very low concentrations, in the millimolar regime.
The least polar compound among them, nitrobenzene, shows the most substantial amplification of binding to the cross-linker at higher concentrations.NB is known to form NB-NB pairs and stacks of the aromatic rings, 69 resulting thus in positive cooperativity for local adsorption (refer to earlier work 38 for further explanation).Note that the bulk concentration of ρ 0 = (28 ± 2) mM might have exceeded the solubility of nitrobenzene in water.At 298 K the experimental value is 16 mM, but computer simulations can overestimate the solubility (115 mM). 70P 0 also performs stacking (Fig. 4c), but due to its additional OH-tail, it has a higher water solubility and is thus less probable to aggregate.In stark contrast, the NP − adsorption to the whole network unit drastically drops at higher concentrations due to their electrostatic repulsion.This is an example of strong negative cooperativity of adsorption at higher concentrations.Note that a real hydrogel may change in size (in particular close to its VPTT) because of the solute-polymer interactions and that solutes may occupy a non-negligible volume, which in return limit the solute adsorption. 71,72,20,35

Partitioning in a hydrogel
The adsorption coefficients retrieved in this work can be used to estimate the resulting solute partition ratios in swollen PNIPAM-BIS hydrogels, which will be compared with experimental data in this section.To this end, we extrapolate our results using an idealized, mean-field model of a large network.

Idealized hydrogels
As a start, the solute partition ratio is determined by the outside bulk concentration ρ 0 and the concentration inside the hydrogel ρ in , namely Here, V gel is the volume of the entire hydrogel including the water and should not be confused with the excluded volume V ex .
Figure 6 Simulation snapshot of NP − (illuminated yellow) benefiting from several possible interaction sites in the BIS-linker (illuminated green) proximity, serving as an illustrative explanation for the strong adsorption amplification due to the cross-linker (see Table 1 and Fig. 5).Nitro-oxygens of NP − can form hydrogen bonds (dashed orange lines depict potential hydrogen bond formation in this configuration) with numerous amide hydrogens, whereas the non-polar aromatic ring is surrounded by hydrophobic environment, i.e., isopropyl groups, (highlighted by bubbles) of the flexible PNIPAM side-chains.An aromatic ring-backbone contact has been observed but less frequently than the presented scenario.NP − can stay in such a conformation (with interchanging binding sites) for several tens of nanoseconds.
The number of particles N in inside the gel can be assessed using the total solute adsorption as If the hydrogel has no net effect, i.e., Γ tot = 0, the concentration inside the gel is equal to the bulk value ρ 0 , and the particles inside the gel account to ρ 0 V gel .Neglecting the effects of dangling (terminal) ends, the adsorption can be assumed to be the sum of single chain adsorption and the effect of all cross-linkers, yielding where N mer and N xlink stand for the number of PNIPAM monomers and BIS-linkers, respectively.Plugging all ingredients into eqn (10), the solute partition ratio can be expressed as with the PNIPAM monomer concentration ρ mer = N mer /V gel , the BIS-to-PNIPAM monomer ratio α = N xlink /N mer , and adsorption coefficients Γ * = Γ/ρ 0 .Considering now a defect-free diamond lattice network architecture 22,49 of the hydrogel, we can deduce the functional form of the monomer concentration vs. the cross-linker ratio, i.e., ρ mer → ρ mer (α), see Appendix B. For different adsorption coefficient pairs (Γ * xlink , Γ * mer ) and in dependence on the cross-linker ratio, K is visualized in Fig. 7.We find that K and α can have a non-linear relation and even a non-monotonic behavior.The reason is that higher crosslinker ratios directly enhance the influence of Γ * xlink , and additionally, as already discussed, increase the PNIPAM concentration ρ mer (α), promoting the influence of Γ * mer .If now Γ * xlink and Γ * mer have even different signs, i.e., a solute, for example, is preferentially desorbed from the polymer but adsorbed by cross-linker then naturally non-monotonic behavior must occur.
Typical values for cross-linker ratios in experiments range from roughly 0.02 to 0.2.In our model (see appendix B), this corresponds to volume fractions ranging from approximately 0.03 to 0.75 with an almost linear relation to α in this interval.Thus the plotted region in Fig. 7 is quite reasonable for demonstrating the non-linear and non-monotonous α-dependencies of the partition ratio.In particular, with positive adsorption coefficients, like all nitro-aromatics have, we find that the partition ratio monotonically increases with larger cross-linker ratios, exemplified by nitrobenzene in Fig. 7.
Selective solute-cross-linker binding affinities are not solely responsible for an increasing partition ratio when increasing the cross-link ratio.As an illustration, we show scenarios of hypothetical solutes that have either zero adsorption to cross-linkers or zero adsorption to the chains.For small values of α, the coefficient Γ * mer has greater impact on partitioning increase than Γ * xlink .In the case of benzene (similarly hexane and butane), where we find a positive chain adsorption, but repulsion from the cross-linker, the partitioning reaches a plateau at α = 0.2.Weaker chain adsorption, or stronger cross-linker repulsion can lead to a maximum in the plotted range, which is exemplified by the hypothetical solute with Γ * mer = 1 nm 3 and Γ * mer = −4 nm 3 .The very opposite case, i.e., cross-linker affinity in combination with chain avoidance, as we find for the tested ion pair Na + /Cl − , exhibits a minimum.
Summing up, the adsorption coefficients Γ * mer and Γ * xlink determine the gradient and concavity of the solute partitioning in dependence on the cross-linker ratio, assuming homogeneous and diamond lattice-like network structure.We conclude that partitioning vs. the cross-linker ratio may be complex and non-monotonous, exhibiting minima and maxima and intercepting the K = 1 line.

Relating to experiments
In real hydrogels, one has to be aware of additional effects.One strong assumption in our model is a rigid and homogeneous network structure, which in general is not the case in the real world.Though techniques have emerged to control the cross-linker density throughout the gel, 73,14,74 the hydrogel structure is still subject to the randomness of the polymerization process and thus retains inhomogeneities.This may lead to nano/micro cavities within the gel and more Dashed lines have the same cross-linker adsorption, solid lines have identical monomer adsorption.Nitrobenzene (NB) has an overall positive binding affinity and K is strictly an increasing function of α.Benzene (B) has a positive adsorption to the chain monomers but a slightly negative cross-linker effect, resulting in a maximum value at α = 0.2.The Na + /Cl − pair is the opposite case, it has a negative chain adsorption coefficient but a positive one for the linker, exhibiting a partitioning minimum.The orange line presents the same case but with zero crosslinker effect (Γ * xlink = 0) and it has a linear relation to the volume fraction.The black dashed line presents the polymer volume fraction φp for this idealized diamond network and its scale is on the right.complex network architectures than assumed, influencing polymer volume fraction and partitioning.
Our investigation focuses on very low solute concentrations, though the response of the polymer to the penetrants might not be negligible.From experiments 71,72,75 it is known, that solutes may change the hydrogel's VPTT, which has additionally been demonstrated in computer simulations. 20,35evertheless, our idealized approach tackling the partitioning does allow for an indirect comparison to experimental data.One experiment on the rate of the nitrobenzene reduction in an (N-isopropylacrylamide-co-acrylic-acid (PNIPAM-co-AAc)) nanoreactor 76 shows an increase in the reduction rate with increasing cross-link (BIS) density, which is attributed to the higher nitrobenzene concentration inside the hydrogel.Parasuraman et al. 77 used a very similar hydrogel (PNIPAM-co-AAc-BIS) and proved the increased dye uptake (Orange II) with increasing cross-link ratio.Both studies qualitatively support our findings.Theym, however, have in common that the initial increase of the reduction rate and the dye uptake, respectively, apparently saturates for higher cross-link degrees.This effect is not captured by our model and might result from steric hindrances, i.e., undersized pore/mesh size of the polymer architecture and already occupied adsorption sites for higher solute concentrations.
After qualitatively confirming the impact of the cross-link ratio on the partitioning of the aromatic compounds, we will now assess the comparison in terms of absolute values.Experimentally, partition ratios have been reported for several molecules containing aromatic rings.A study by Molina et al. 52 retrieved K in PNIPAM-BIS-hydrogels (α = 2%) for probe drugs (tryptophan, propranolol chloride, dansyl chloride, methyl orange, riboflavin and ruthenium-tris(2,2'bibyridiyl) dichloride), which contain two to six aromatic rings as well as polar and/or charged residues.The partition ratio, depending on the compound, ranges from roughly 4.6 to 10.
Comparing to our perfect network model, the much smaller NB and NP − show partition ratios of approximately 1.3 and 1.8 respectively at α = 2%, and 2.7 and 5.7 at α = 5% in the low dilution limit (N solute = 1).For N solute = 20, K is about 3.7 at α = 2% for NB.It is expected, that larger molecules at higher concentrations, as established in the mentioned experiments, will lead to higher partition ratios 78 and can thus be regarded as supportive of our results.Furthermore, the adsorption increase due to positive cooperativity (e.g., NB) has been shown for methylene blue in a superabsorbent hydrogel. 79he salt partition ratios in 1% cross-linked PNIPAM-BIS gels at room temperature have been reported 80 and amount to K LiCl = 0.97 ± 0.05, K KCl = 0.91 ± 0.05 and K NaCl ≈ 0.95, i.e., just below unity as in our prediction for sodium chloride.
However, the non-monotonicity of the partition ratio in dependence on the cross-link ratio predicted by our model has not been reported by experimentalists so far and is yet to be tested.

Concluding remarks
We investigated the effects of cross-linking on solute adsorption in swollen PNIPAM hydrogels by means of explicitwater MD simulations at T = 290 K, i.e., below the PNI-PAM collapse transition temperature.We considered a generic hydrogel subunit consisting of one BIS-linker and four PNIPAM chains, which was kept in tetrahedral geometry by position-restrained backbone terminals.By subdividing the radial distance from the central cross-linker, we classified four different adsorption regions according to the polymer's prevalent features, namely the cross-linker region, the linear chain region, the chain terminal region, and the bulk solvent domain.We evaluated the adsorption of different solutes, representing typical charged, polar and nonpolar molecular compounds.
Comparing the cross-linker and monomer effects on adsorption, we find different scenarios.Apolar species show small attraction to chain monomers and slight repulsion from the cross-linker region.Sodium chloride behaves the opposite, it has negative chain adsorption but is attracted towards the cross-linker.The strongest adsorbing solutes, the nitro-aromatics, adsorb to all parts of the polymer and show the highest binding affinity, which is promoted by hydrophobic interactions between the aromatic ring and PNI-PAM's isopropyl groups as well as by hydrogen bonds be-tween the nitro-oxygens and the amide groups.The adsorption at the cross-linker relative to a single PNIPAM monomer adsorption, spans from Γ xlink /Γ mer ≈ 2 (NB) to Γ xlink /Γ mer ≈ 10 (NP 0 , NP − ).This indicates that the crosslinker can significantly enhance the overall adsorption to the network unit.Hence, for solutes that have a significant affinity to PNIPAM chains already, the dense cross-linker region, where many-body attractions are at play, amplifies the local adsorption by even an order of magnitude.Thereby we confirm the 'vertex trapping' effect that has been first reported in generic coarse-grained simulations of polymer networks. 24,25,26,27,28n the case of the nitro-aromatics, we furthermore performed simulations with higher solute concentrations to estimate cooperative adsorption effects.Nitrobenzene shows enhanced aromatic stacking at the cross-linker and the adsorption is elevated in a superlinear fashion with increasing concentration.Nitrophenol shows similar, positive cooperativity but a less pronounced behavior.For both NB and NP 0 the cross-linker promotes higher positive cooperation effects than single chains.In contrast, in the case of the negatively charged nitrophenolate, we observe less adsorption at higher concentrations due to negative cooperativity stemming from the electrostatic repulsion.
The adsorption coefficients for cross-linker and chain monomers in the low concentration regime were used to estimate partition ratios of the solutes within an idealized, homogeneous diamond-lattice macrogel, which allowed a comparison with experimental findings.In our model, we found that highly adsorbing substances like nitro-aromatics have a partition ratio ranging from 2 to 5 at a cross-linker concentration of 5%.Solutes with adsorption coefficients of opposite signs may show non-monotonic behavior as a function of the cross-linker ratio: Positive/negative chain adsorption and negative/positive cross-linker adsorption leads to a maximum/minimum in the partition ratio.These yet poorly known features should be considered in future experiments and modeling of hydrogels as they play an important role for the fine-tuning of solute uptake within the needs of the desired function and application.

A Force field parameters
From our ab-initio calculations with the Gaussian 09 software, 47 we obtain partial charges for the cross-linker and monomers, which are very similar and in good agreement with the OPLS-QM2 force field. 31Due to the chemical similarity of BIS and PNIPAM, we apply the PNIPAM's partial charges of the OPLS-QM2 force field for BIS atoms as well.The partial charges from the original OPLS-AA 46 and OPLS-QM2 force fields as well as our results are shown in Table 2.The standard OPLS force field is used for all remaining parameters like potentials, masses, etc.The hitherto undefined N-C-N angle and C-N-C-N dihedral potentials for the cross-linker were adopted from the OPLS C-C-C angle and C-N-C-C dihedral parameters, respectively.Recapturing the results for the adsorptions around single chains reported by our group before, 38 the reader may compare the standard OPLS-AA force field for PNIPAM polymers with the further optimized OPLS-QM2 version 31 employed in this work.The qualitative trends, i.e., size and polarity dependence as well as the strong nitro aromatic binding affinity remain.However, comparing absolute numbers for aromatic compounds, the standard force field 38 shows a roughly twice as high adsorption coefficient to the chain monomers.We attribute this effect to more polar amide potentials, masses, etc.The so far undefined N-C-N-angle and N-C-N-C-dihedral potentials for the crosslinker were adopted from the OPLS C-C-C-angle and C-N-C-C-dihedral parameters.  2 .
Eventually, the monomer concentration as well as the solute partition coefficient, eqn (13), are fully defined.By knowing n mer (α) we can further approximate the polymer volume fraction, which allows us to identify physically meaningful values of α and K. Using the water profile (Fig. 3) we extract the excluded volume per monomer as V ex,mer = 0.167 nm 3 and assume the BIS' volume to scale with the number of heavy atoms (without hydrogens) compared to a PNIPAM monomer (11. vs. 8), yielding V ex,xlink = 0.230 nm 3 .The excluded volume in one unit cell is then v ex = 8V ex,xlink + n mer (α)V ex,mer and the  2.
groups in the OPLS-QM2 force field and hence to more hydrophilic behavior.Moreover, the adsorption of NP − to PNIPAM (in the OPLS-QM2 version) is significantly larger owing to its polar nitro group.

B Defect-free PNIPAM-BIS hydrogels
We assume a homogeneous network with cross-linkers arranged in a perfect diamond lattice and thus connected with chains of an equal number of monomers.Given this perfect lattice, we can simply deduce from the geometry the monomer concentration as a function of the cross-linker ratio, ρ mer (α) = n mer (α)/v gel (α).We can also assess the network's unit cell volume v gel (α) and its PNIPAM monomer content n mer (α).Such a unit cell contains eight crosslinkers and 16 chains, hence n mer (α) = 8/α, and the number of monomers between two associated cross-linkers is n chain (α) = 1/(2α).The distance between two cross-linkers is then given by xx (α) = xlink + n chain (α)λ∆L c , (14)   where xlink = 0.22 nm, the effective length contribution of one BIS-linker, which is an average value retrieved from simulations solving L = xlink /2+ ee .The unit cell volume reads v gel (α) = 8 with θ ≈ 109.5 • being the angle between any pair of adjacent chains.Eventually, the monomer concentration as well as the solute partition ratio, eqn (13), are fully defined.By knowing n mer (α), we can further approximate the polymer volume fraction, which allows us to identify physically meaningful values of α and K. Using the water profile (Fig. 3) we extract the excluded volume per monomer as V ex,mer = 0.167 nm 3 and assume the BIS' volume to scale with the number of heavy atoms (without hydrogens) compared to a PNIPAM monomer (11. vs. 8), yielding V ex,xlink = 0.230 nm 3 .The excluded volume in one unit cell is then v ex = 8V ex,xlink + n mer (α)V ex,mer and the polymer volume fraction for different values of α is obtained as φ p = v ex /v gel (Fig. 7).

Fig. 1
Fig. 1 Chemical structures (top) and corresponding ball-and-stick representations (bottom) of poly-(N-isopropylacrylamide) [PNIPAM] (panels a, b), and the crosslinker N,N -methylenebisacrylamide [BIS] (panels c, d).Dashed lines represent bonds to neighboring PNIPAM monomers.Associated carbon atoms are referred to as the polymer backbone.PNIPAM's amide group and isopropyl group are called side chains.Amide groups and BIS central carbon atoms are hydrophilic, potentially forming hydrogen bonds with the surrounding.The backbone and isopropyl group have a hydrophobic character.

Figure 1
Figure 1 Chemical structures (top) and corresponding ball-and-stick representations (bottom) of poly(N -isopropylacrylamide) [PNIPAM] (panels a, b), and the cross-linker N ,N -methylenebisacrylamide [BIS] (panels c, d).Dashed lines represent possible bonds to neighboring PNIPAM monomers or BIS.Associated carbon atoms are referred to as the polymer backbone.PNIPAM's amide group and isopropyl group form side chains.The amide groups and BIS' central methylene bridge are hydrophilic, potentially forming hydrogen bonds 45 with the surrounding.The backbone and isopropyl group have a hydrophobic character.

Figure 2
Figure 2 Simulation snapshots of the studied polymeric molecule, consisting of one BIS-linker and four PNIPAM chains, at T = 290 K and relative chain extension λ = 0.83 in (a) all-atom ball-and-stick representation and (b) showing only backbone and heavy BIS atoms inside the simulations box (--) and in addition as 2D projections.Terminal backbone carbons are position-restrained (marked as ×) in the corners of a virtual tetrahedron (---), which has an edge length of 4.69 nm and a vertex-to-centroid distance of L = 2.87 nm.For the sake of clarity, water molecules are not shown.(c) Probe solutes.Atoms are color-coded by element, i.e., black carbon, red oxygen, blue nitrogen, and white hydrogen atoms.See Table1for the full simulation specifications.

Figure 3
Figure 3 The different adsorption domains I, II, and III resolved in the radial distance r to the COM of the cross-linker, illustrated in (a) a simulation snapshot of BIS plus two PNIPAM branches and (b) by radial water density (left hand axis) and PNIPAM monomer number (right hand axis) profiles.The first interval I is the 'cross-link region', where BIS as well as collective PNIPAM effects mix.The PNIPAM profile shows shortrange oscillations, leading to the non-monotonicity in the water profile.The range of influence of the cross-linker is assumed to vanish around ra = 1.5 nm.Region II, ranging from ra to r b , is dominated by the linear PNIPAM chain.Here, the monomer number, ρmer(r)4πr 2 , is roughly constant.In interval III, from r b = 2.5 nm to rc = 3.5 nm, the polymer chains terminate.The average distance between the position-restrained backbone terminals and the cross-linker, L = 2.87 nm, is shown by the vertical dotted (red) line.For large distances, i.e., r > rc, we assume negligible influence from the polymer and consider region IV as 'bulk'.

Figure 4
Figure 4 (a) Normalized radial density profiles of the solutes with respect to the COM of the BIS-linker for the low-adsorption species (C 6 , C 4 , C 1 OH, Na + , Cl − ), plus benzene (--B) and water (--H 2 O) profiles for comparison.Note that in the case of methanol (--C 1 OH) and sodium chloride (• • • • Na + /Cl − ) simulations were carried out with 30 molecules or 3 ion pairs, respectively.Panel (b) shows much stronger adsorbing molecules BZ, NB, NP 0 , NP − , all of aromatic nature, in the high dilution limit (--) and more concentrated solutions with a total of 20 (---) molecules (and 20 Na + counterions in the case of NP − ) per simulation.(c) Simulation snapshot showing a stacking of NP 0 in the cross-link region.BIS and NP 0 are highlighted green and yellow respectively, PNIPAM is shown in licorice representation.Note that the aggregation of the solutes usually looks less ordered than presented due to thermal fluctuations.

Figure 5
Figure 5 Adsorption coefficients (Γ * = Γ/ρ 0 ) quantifying the binding affinities of different solutes to one PNIPAM monomer Γ * mer and the crosslinker Γ *xlink in the infinite dilution regime, as summarized in Table1.The adsorption to the monomers increases with size, e.g., compare C 4 → C 6 .The aromatic molecules show the highest adsorption, which is magnified by increasing polarity B → NB → NP 0 → NP − , coinciding with a crosslinker-enhanced binding affinity.

Figure 7
Figure 7Partition ratios for different solutes in an ideal diamond-lattice polymer network (see appendix B) as a function of cross-linker ratio.The curves result from a competition between the adsorption coefficients Γ * mer and Γ * xlink .Dashed lines have the same cross-linker adsorption, solid lines have identical monomer adsorption.Nitrobenzene (NB) has an overall positive binding affinity and K is strictly an increasing function of α.Benzene (B) has a positive adsorption to the chain monomers but a slightly negative cross-linker effect, resulting in a maximum value at α = 0.2.The Na + /Cl − pair is the opposite case, it has a negative chain adsorption coefficient but a positive one for the linker, exhibiting a partitioning minimum.The orange line presents the same case but with zero crosslinker effect (Γ * xlink = 0) and it has a linear relation to the volume fraction.The black dashed line presents the polymer volume fraction φp for this idealized diamond network and its scale is on the right.

C 2 H 2 C 1 3 C 5 H 5 3 C 2 H 2 C 1 Fig. 8
Fig. 8 Chemical structures of PNIPAM and BIS.The superscript indices differentiate atoms regarding their partial charges, shown in Table2.

9 Figure 8
Figure 8 Chemical structures of PNIPAM and BIS.superscript indices differentiate atoms regarding their partial charges, shown in Table2.

Table 2
Partial charges (in unit charges) of different force fields and results from own ab-initio (Hartree-Fock) calculations on the 6-31G(d)-level with electrostatic potential fitting for BIS and PNIPAM.Chemical structure with the different atom classification are depicted in Fig. 8. BIS charges of our calculations were retrieved by analyzing 158 conformations of a molecule consisting of one BIS cross-linker and four PNIPAM monomers attached.PNIPAM charges were obtained from 36 different single chains consisting of 11 monomers.In this work we use the OPLS-QM2 partial charges for PNIPAM and BIS.The charges for C 6 and H 6 atoms in BIS were chosen similar to the C 4 and H 4 charges in PNIPAM and with overall electroneutrality of the cross-linker in mind.A small test simulation with slightly redistributed partial charges of the central C 6 H 6 2 -group (C 6 = 0.1,H 6 = 0.15,C 3 = 0.56) did not change the nitrobenzene adsorption affinity.The overall positive partial charges of these atoms (C 6 H 6