Molecular Dynamics Simulation of Hydrogels Based on Phosphorylcholine-Containing Copolymers for Soft Contact Lens Applications

The structure and dynamics of copolymers of 2-hydroxyethyl methacrylate (HEMA) with 2-methacryloyloxyethyl phosphorylcholine (MPC) were studied by molecular dynamics simulations. In total, 20 systems were analyzed. They differed in numerical fractions of the MPC in the copolymer chain, equal to 0.26 and 0.74, in the sequence of mers, block and random, and the water content, from 0 to 60% by mass. HEMA side chains proved relatively rigid and stable in all considered configurations. MPC side chains, in contrast, were mobile and flexible. Water substantially influenced their dynamics. The copolymer swelling caused by water resulted in diffusion channels, pronounced in highly hydrated systems. Water in the hydrates existed in two states: those that bond to the polymer chain and the free one; the latter was similar to bulk water but with a lower self-diffusion coefficient. The results proved that molecular dynamics simulations could facilitate the preliminary selection of the polymer materials for specific purposes before their synthesis.

Hydroxyl and carbonyl functional groups in HEMA molecules (Figure 1) can hydrogenbond water molecules.On the other hand, hydrophobic methyl groups in the polymer matrix ensure hydrolytic stability and mechanical strength.Moreover, P-HEMA is inert to many chemicals, resistant to degradation, thermally stable, and has good mechanical properties as an SCL material due to cross-linking of the chains [11][12][13][14][15].The main disadvantage of pure P-HEMA hydrogels in ophthalmologic applications is low equilibrium water content (EWC) of only 38% by mass.This results in poor oxygen permeability as oxygen is transported mainly through the water phase of the hydrogels [16].Consequently, long-term hypoxia and hypercapnia (increased level of carbon dioxide) may occur with the use of P-HEMA SCLs and lead to deterioration of health [17].Copolymers of HEMA with compounds that bind more water are better in this respect [18][19][20].
High water-capacity units in the copolymer can be those of 2-methacryloyloxyethyl phosphorylcholine, MPC (Figure 1).The phosphorylcholine (PC) moiety enhances bioinertia as a consequence of impaired electrostatic interactions and the formation of a hydration shell with a large amount of highly mobile free water around [21][22][23][24][25][26].The latter strongly repulses the surrounding proteins [27,28].Moreover, the PC moiety inhibits cell adhesion and reduces blood coagulation, and improves immune response [27,29,30].These properties make MPC an ideal component of SCL biomaterials [31][32][33][34][35].However, the properties of copolymers consisting of the same constituent units may vary depending on their arrangement in the chain.The chain arrangement may be complex even for linear biopolymers, which could be classified as statistical, random, or alternating.In the statistical arrangement, the sequential distribution obeys statistical laws, e.g., Markovian statistics, while in the random distribution, the probability of a given mer occupying a site in the chain does not depend on the adjacent mers.Block copolymers are composed of block macromolecules built of many constitutional units each [36].High water-capacity units in the copolymer can be those of 2-methacryloyloxyethyl phosphorylcholine, MPC (Figure 1).The phosphorylcholine (PC) moiety enhances bioinertia as a consequence of impaired electrostatic interactions and the formation of a hydration shell with a large amount of highly mobile free water around [21][22][23][24][25][26].The latter strongly repulses the surrounding proteins [27,28].Moreover, the PC moiety inhibits cell adhesion and reduces blood coagulation, and improves immune response [27, 29,30].These properties make MPC an ideal component of SCL biomaterials [31][32][33][34][35].However, the properties of copolymers consisting of the same constituent units may vary depending on their arrangement in the chain.The chain arrangement may be complex even for linear biopolymers, which could be classified as statistical, random, or alternating.In the statistical arrangement, the sequential distribution obeys statistical laws, e.g., Markovian statistics, while in the random distribution, the probability of a given mer occupying a site in the chain does not depend on the adjacent mers.Block copolymers are composed of block macromolecules built of many constitutional units each [36].
The polymer chain arrangement has a significant influence on the physicochemical properties of the materials but it is difficult to test experimentally.On the other side, despite significant progress in improving the properties of the materials used in the production of SCLs, searching for new biomaterials for advanced ophthalmological applications is still a matter of interest.Computer simulations, complementary to experimental studies, predicting the physicochemical properties of polymers may substantially shorten the time necessary for designing new polymer materials and give a new look at their properties.Nevertheless, reports on computer simulations of the physicochemical properties of polymers potentially applicable in ophthalmology are The polymer chain arrangement has a significant influence on the physicochemical properties of the materials but it is difficult to test experimentally.On the other side, despite significant progress in improving the properties of the materials used in the production of SCLs, searching for new biomaterials for advanced ophthalmological applications is still a matter of interest.Computer simulations, complementary to experimental studies, predicting the physicochemical properties of polymers may substantially shorten the time necessary for designing new polymer materials and give a new look at their properties.Nevertheless, reports on computer simulations of the physicochemical properties of polymers potentially applicable in ophthalmology are scarce [37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52].Therefore, the present work was aimed at filling this gap at least partially.The attention was focused on the structural and dynamic properties of the P(MPC-co-HEMA) hydrogels, especially on the hydration phenomena in these materials.We studied the mechanism of water diffusion which influences the transport of oxygen through the hydrogels.Proper oxygen concentration on the eyeball surface is crucial for the health and comfort of the SCL wearers.Understanding transport processes in polymeric materials facilitates an assessment of their suitability for ophthalmological applications.

Results and Discussion
In the simulations, we considered several intermolecular distances and interactions.To keep the text readable, we used shortcuts rather than full descriptions of the atoms and interactions.Thus, we marked the atoms by the subscripts H, M, and W, which denote HEMA, MPC, and water molecules, respectively.The full identification of the symbols for HEMA and MPC is given in Figure 1.
Copolymers of 2-hydroxyethyl methacrylate (HEMA) with 2-methacryloyloxyethyl phosphorylcholine (MPC), P(MPC-co-HEMA), dry and hydrated with water content 10, 20, 40, and 60% by mass were investigated.The P(MPC-co-HEMA) copolymer chains in the simulations were composed of 50 mers, 13 of one kind and 37 of the other, in block or random configuration.The following four copolymers were studied: random with 13 MPC and 37 HEMA mers, HMH6MH4MH4MH4MHMHMHMH11MHM2H2M2H (R13), • random with 37 MPC and 13 HEMA mers, M10HM4H2M6H3M5H4MHM5HM3HM3 (R37), where M stands for MPC and H for HEMA, the subscripts denote the number of mers, and codes in parentheses identify the copolymer.From now on, the codes will be used for simplicity.

Structure and Chain Packing
The flexibility of the HEMA and MPC mers would be manifested in differences in their geometry resulting from the copolymer chain arrangement and degree of hydration.We chose distances between two atoms in each of the mers as indicators of this property.For both mers, the first atom was that of carbon belonging to the copolymer main chain, i.e., C H or C M in Figure 1.The second atom was located at the end of the side chain: the O H3 of HEMA and N M of MPC.
First, the HEMA chains will be discussed.The radial distribution functions g(r) of the C H -O H3 pairs are shown in Figure 2. The shortest distance between C H and O H3 atoms assessed from the RDF first peak maximum is 5.0 Å for all the systems studied.Thus, neither the arrangement of the copolymer chain nor the hydration affects the HEMA geometry.This proves that the HEMA mer is relatively rigid.Higher g(r) peaks for hydrated polymers result from the lower number density of HEMA mers in cells enlarged for accommodating water molecules (cf.Table 1).The MPC groups are different in this respect.Figure 3 shows the RDFs calculated for the CM-NM pair.For non-hydrated copolymers, the CM-NM interatomic distance is equal to 5.7 Å for R13 and 5.6 Å for B13.It is significantly shorter for the copolymers with 37 mers of MPC in the chain: 4.9 Å and 4.3 Å for R37 and B37, respectively.Thus, the MPC chains are stiffened by those of HEMA in their vicinity, while they can coil when in bigger blocks formed by MPC mers.The CM-NM distances shorter than the CH-OH3 ones seem rather odd because the MPC mer is almost twice as long as that of HEMA.An explanation could be the coulombic attraction between oppositely charged trimethylammonium and phosphate groups of adjacent mers causing twisting and shortening of the side chain of the polymer (Figure 4).The MPC groups are different in this respect.Figure 3 shows the RDFs calculated for the C M -N M pair.For non-hydrated copolymers, the C M -N M interatomic distance is equal to 5.7 Å for R13 and 5.6 Å for B13.It is significantly shorter for the copolymers with 37 mers of MPC in the chain: 4.9 Å and 4.3 Å for R37 and B37, respectively.Thus, the MPC chains are stiffened by those of HEMA in their vicinity, while they can coil when in bigger blocks formed by MPC mers.The C M -N M distances shorter than the C H -O H3 ones seem rather odd because the MPC mer is almost twice as long as that of HEMA.An explanation could be the coulombic attraction between oppositely charged trimethylammonium and phosphate groups of adjacent mers causing twisting and shortening of the side chain of the polymer (Figure 4).
The C M -N M distance increases due to the hydration of the polymer, for R37 and B37 in particular, and eventually reaches 6.3 Å for the four systems with 60% of water.However, the MPC chains remain curved.For comparison, the distance in a fully straightened chain would be about 11.2 Å, as estimated from the molecule geometry and bond lengths.The elongation of C M -N M distance suggests that water penetrates the voids between the MPC side chains and influences their shape.Two mechanisms could lead to this effect.First, non-specific hydration of the trimethylammonium and phosphate groups may weaken the coulombic attraction between the adjacent chains because of the high dielectric constant of water.Second, water bridges hydrogen-bonded to the chains may stabilize their geometry.This will be discussed in the next sections.
Although some structural features are manifested in the reported RDFs, the lack of long-range regularity evidences the amorphousness of the studied polymers.The CM-NM distance increases due to the hydration of the polymer, for R37 and B37 in particular, and eventually reaches 6.3 Å for the four systems with 60% of water.The CM-NM distance increases due to the hydration of the polymer, for R37 and B37 in particular, and eventually reaches 6.3 Å for the four systems with 60% of water.

Distribution of Water Molecules
In this section, the location of water molecules in the hydrogels will be discussed in detail beginning from the vicinity of HEMA chains.The RDF functions of the C H -O W , O H3 -O W , and O H2 -O W pairs calculated for the hydrated R37 copolymer are plotted in Figure 5.
bond lengths.The elongation of CM-NM distance suggests that water penetrates the voids between the MPC side chains and influences their shape.Two mechanisms could lead to this effect.First, non-specific hydration of the trimethylammonium and phosphate groups may weaken the coulombic attraction between the adjacent chains because of the high dielectric constant of water.Second, water bridges hydrogen-bonded to the chains may stabilize their geometry.This will be discussed in the next sections.
Although some structural features are manifested in the reported RDFs, the lack of long-range regularity evidences the amorphousness of the studied polymers.

Distribution of Water Molecules
In this section, the location of water molecules in the hydrogels will be discussed in detail beginning from the vicinity of HEMA chains.The RDF functions of the CH-OW, OH3-OW, and OH2-OW pairs calculated for the hydrated R37 copolymer are plotted in Figure 5.The RDF functions for B13, B37, and R13 are similar to the presented ones.The shortest distance between OH3 and OW is 3.1 Å and it does not depend on the amount of water in the system, at least within the concentration range studied.Similarly, a short distance equal to 3.3 Å is between OH2 and OW, but the respective g(r) peak is tiny, while a pronounced maximum occurs at 5.2 Å.The two ca. 3 Å-long distances suggest hydrogen bonds occurring between water molecules and the ether or terminal hydroxyl groups of the HEMA chain.The second interaction predominates over the first most probably because of the main polymer chain hindering the access of the water OH group to the carbonyl oxygen atom.The OH2-OW distance of 5.2 Å, corresponding to the pronounced The RDF functions for B13, B37, and R13 are similar to the presented ones.The shortest distance between O H3 and O W is 3.1 Å and it does not depend on the amount of water in the system, at least within the concentration range studied.Similarly, a short distance equal to 3.3 Å is between O H2 and O W , but the respective g(r) peak is tiny, while a pronounced maximum occurs at 5.2 Å.The two ca. 3 Å-long distances suggest hydrogen bonds occurring between water molecules and the ether or terminal hydroxyl groups of the HEMA chain.The second interaction predominates over the first most probably because of the main polymer chain hindering the access of the water OH group to the carbonyl oxygen atom.The O H2 -O W distance of 5.2 Å, corresponding to the pronounced peak at g(r), and that of 5.0 Å between C H and O W are typical of non-hydrogen-bonded molecules.The shortest distance between PM and OW, 4.1 Å, suggests that water molecules gather in the vicinity of the negatively charged phosphate group.Similar gathering around the positive trimethylammonium group also occurs, but at a longer distance of 4.7 Å because of the three blocking methyl groups.The distance between CM and OW, 5.0 Å, is equal to that between CH and OW.Indeed, the CH and CM atoms belong to the main polymer chain, and the HEMA and MPC side chains are identical up to the second O atom (OH3 and OM3).

Hydrogen Bonding
To evaluate the feasibility of water-polymer hydrogen bonds, we considered all potential hydrogen-bonding sites: the oxygen, nitrogen, and phosphorus atoms in HEMA and MPC mers, as well as hydrogen atoms from water molecules.As in previous sections, we will begin the discussion from the HEMA mers.
The RDF functions for the three pairs (OH1 or OH2 or OH3 in R37)-HW are shown in Figure 7.The distances OH1-HW and OH3-HW corresponding to the RDFs' first maxima are equal to 2.2 Å suggesting hydrogen bonds.However, the hydroxyl groups (OH3) are probably donors of the hydrogen atoms rather than acceptors.That could be concluded because the first peak of the respective RDF is just a small local maximum preceding a steep increase up to the global maximum at r ≈ 3 Å for the systems with water content not exceeding 40% by mass.The OH3-H…OW bonds seem very probable because of the short OH3-OW distance equal to 3.1 Å (see Figure 8).Contrary to OH3, the ether and carbonyl oxygen atoms, OH2 and OH1, can only be acceptors in hydrogen bonding.From these two, OH1 is doubtless the stronger acceptor that is manifested in prominent first maxima at r = 3.1 Å and the long-range oscillations of the OH1-OW RDFs.The oscillations would be

Hydrogen Bonding
To evaluate the feasibility of water-polymer hydrogen bonds, we considered all potential hydrogen-bonding sites: the oxygen, nitrogen, and phosphorus atoms in HEMA and MPC mers, as well as hydrogen atoms from water molecules.As in previous sections, we will begin the discussion from the HEMA mers.
The RDF functions for the three pairs (O H1 or O H2 or O H3 in R37)-H W are shown in Figure 7.The distances O H1 -H W and O H3 -H W corresponding to the RDFs' first maxima are equal to 2.2 Å suggesting hydrogen bonds.However, the hydroxyl groups (O H3 ) are probably donors of the hydrogen atoms rather than acceptors.That could be concluded because the first peak of the respective RDF is just a small local maximum preceding a steep increase up to the global maximum at r ≈ 3 Å for the systems with water content not exceeding 40% by mass.The O H3 -H. ..O W bonds seem very probable because of the short O H3 -O W distance equal to 3.1 Å (see Figure 8).Contrary to O H3 , the ether and carbonyl oxygen atoms, O H2 and O H1 , can only be acceptors in hydrogen bonding.From these two, O H1 is doubtless the stronger acceptor that is manifested in prominent first maxima at r = 3.1 Å and the long-range oscillations of the O H1 -O W RDFs.The oscillations would be reinforced by the cooperativity of hydrogen bonds.However, the cooperativity could not be detected because the electron structures of the modeled molecules remained unchanged during the simulation.Of course, increased water concentration causes more water molecules to stay away from the oxygen atoms of HEMA and the g(r) maxima to decrease.That has been reported for other hydrated polymers as well [53,54].
reinforced by the cooperativity of hydrogen bonds.However, the cooperativity could not be detected because the electron structures of the modeled molecules remained unchanged during the simulation.Of course, increased water concentration causes more water molecules to stay away from the oxygen atoms of HEMA and the g(r) maxima to decrease.That has been reported for other hydrated polymers as well [53,54].MPC is a derivative of HEMA with phosphocholine substituted for hydrogen in the terminal hydroxyl group, cf. Figure 1.Thus, the investigation of the MPC hydration shell was focused on two questions: (i) How are water molecules located around the phosphocholine group, and (ii) Is the arrangement in the remaining area different from that around HEMA?Only the results for R37 are reported.Those for the three other copolymers are similar.The RDFs for six atom pairs, OM1,…, M6-HW, as well as NM-Hw and PM-HW, are plotted in Figures 9 and 10.Respective functions for the pairs with OW rather MPC is a derivative of HEMA with phosphocholine substituted for hydrogen in the terminal hydroxyl group, cf. Figure 1.Thus, the investigation of the MPC hydration shell was focused on two questions: (i) How are water molecules located around the phosphocholine group, and (ii) Is the arrangement in the remaining area different from that around HEMA?Only the results for R37 are reported.Those for the three other copolymers are similar.The RDFs for six atom pairs, O M1,. .., M6 -H W , as well as N M -Hw and P M -H W , are plotted in Figures 9 and 10.Respective functions for the pairs with O W rather than H W are reported in Figures 11 and 12.The RDFs for the carbonyl and ether oxygen atoms in MPC, O M1 , and O M2 , do not differ significantly from the respective functions obtained for HEMA.This suggests that hydration spheres around these fragments of the two sidechains are similar to one another.Sharp peaks on the RDFs, these (see Figure 9) and those for PM-OW (see Figure 12) and PM-HW (see Figure 10), confirm a high affinity of water to the phosphate group.Thus, water forms hydrogen bonds with the OM4 and OM5 atoms that result in the RDF first peak coordinates of 3.1 Å and 2.2 Å for the OM(4,5)-OW and OM(4,5)-HW pairs, respectively.The first peak coordinates for OM3 and OM6-water pairs are equal to 3.1 Å and 2.2 Å too, but the peaks are significantly lower.This suggests the formation of hydrogen bonds involving the OM4 and OM5 atoms facilitated by the negative charge of the phosphate group attracting water molecules through the ion-dipole interaction.We should keep in mind that hydrogen bonding is a purely electrostatic interaction in this simulation.For OM3 and OM6, the effect is weaker if any.
The affinity of water to the nitrogen atom of the choline group, NM, is strong, as the       Sharp peaks on the RDFs, these (see Figure 9) and those for P M -O W (see Figure 12) and P M -H W (see Figure 10), confirm a high affinity of water to the phosphate group.Thus, water forms hydrogen bonds with the O M4 and O M5 atoms that result in the RDF first peak coordinates of 3.1 Å and 2.2 Å for the O M(4,5) -O W and O M(4,5) -H W pairs, respectively.The first peak coordinates for O M3 and O M6 -water pairs are equal to 3.1 Å and 2.2 Å too, but the peaks are significantly lower.This suggests the formation of hydrogen bonds involving the O M4 and O M5 atoms facilitated by the negative charge of the phosphate group attracting water molecules through the ion-dipole interaction.We should keep in mind that hydrogen bonding is a purely electrostatic interaction in this simulation.For O M3 and O M6 , the effect is weaker if any.
The affinity of water to the nitrogen atom of the choline group, N M , is strong, as the first peaks on the RDFs for N M -O W (see Figure 12) and N M -H W (see Figure 10) are high and sharp.However, the hydrogen bonding between N M and water does not occur because of the steric hindrance by the three methyl groups.Thus, the ion-dipole attraction is the main force that keeps water molecules in the vicinity of the group.Doubtless, the preferred orientation of a water molecule is that with the pairs of lone electrons directed towards the positively charged nitrogen atom.It is manifested in the RDFs maxima at 4.7 Å for N M -O W and at 5.0 Å for N M -H W .

Free Volumes
The fractional free volumes (FFV) for twenty P(MPC-co-HEMA) systems determined for a probe with a radius of 0 Å (cf.Section 3.2) are reported in Table 2.The considered polymer material consisted of the copolymer itself and, if it was hydrated, water as well.Figure 13 shows an illustration for R37.Morphologies for the other polymer systems are similar to these, albeit the FFVs for R37 and B37 are slightly bigger than those for R13 and B13 at low water contents in particular.The sequence of mers (block or random) does not affect the FFV significantly.The FFVs decrease with the increase in water content (Table 2), evidencing that swelling does not compensate for the filling of the voids by water.Moreover, the average void shrinks when the water content increases.In real systems, the free volumes for an infinitesimal probe are probably less important than those for molecules of finite size.Figure 14 shows the FFV of non-hydrated P(MPC-co-HEMA) copolymers tested with probes of various radii.The FFV is the smallest in B13 independently of the probe size.For the probes of the "molecular" size, i.e., those of radii above 1.5 Å, the sequence of polymers according to the FFVs is the following: B13 < B37 < R13 ≈ R37.This confirms the intuition-in random arrangements, the phosphocholine groups of long MPC side chains are separated one from another by broader spaces than in the block configurations.That makes penetration of the interchain spaces easier.In real systems, the free volumes for an infinitesimal probe are probably less important than those for molecules of finite size.Figure 14 shows the FFV of non-hydrated P(MPC-co-HEMA) copolymers tested with probes of various radii.The FFV is the smallest in B13 independently of the probe size.For the probes of the "molecular" size, i.e., those of radii above 1.5 Å, the sequence of polymers according to the FFVs is the following: B13 < B37 < R13 ≈ R37.This confirms the intuition-in random arrangements, the phosphocholine groups of long MPC side chains are separated one from another by broader spaces than in the block configurations.That makes penetration of the interchain spaces easier.
Another problem is an influence of swelling on the volume available for water in the polymers.To solve it, the calculations for a hard-sphere probe with a radius of 1.4 Å were performed.In these calculations, the polymers' geometries remained as they were when swelled, while the hydration water was removed.The results are reported in Table 3.In dry copolymers, the FFVs available for water are bigger for random than block configuration.From among the four dry polymers, B13 shows the smallest fractional free volume of up to 20% of water by mass.However, the differences are small and vanish for higher water content.P(MPC-co-HEMA) copolymers tested with probes of various radii.The FFV is the smallest in B13 independently of the probe size.For the probes of the "molecular" size, i.e., those of radii above 1.5 Å, the sequence of polymers according to the FFVs is the following: B13 < B37 < R13 ≈ R37.This confirms the intuition-in random arrangements, the phosphocholine groups of long MPC side chains are separated one from another by broader spaces than in the block configurations.That makes penetration of the interchain spaces easier.Another problem is an influence of swelling on the volume available for water in the polymers.To solve it, the calculations for a hard-sphere probe with a radius of 1.4 Å were performed.In these calculations, the polymers' geometries remained as they were when swelled, while the hydration water was removed.The results are reported in Table 3.In  Examples of the free volume geometry in R37 are drawn in Figure 15.As the polymer hydration increases, the free volumes around the polymer chains become bigger.Eventually, pronounced diffusion channels arise.In the 40 and 60% hydrates, water is not only adsorbed into the polymer matrix but also accumulated around the rolled chains.Such "free" water would not interact directly with the polymer chains but resemble bulk water.This conclusion agrees with the experimental results evidencing bulk water in the hydrates containing at least 30% of water by mass [52,55,56].
Molecules 2023, 28, x FOR PEER REVIEW 14 of 23 dry copolymers, the FFVs available for water are bigger for random than block configuration.From among the four dry polymers, B13 shows the smallest fractional free volume of up to 20% of water by mass.However, the differences are small and vanish for higher water content.Examples of the free volume geometry in R37 are drawn in Figure 15.As the polymer hydration increases, the free volumes around the polymer chains become bigger.Eventually, pronounced diffusion channels arise.In the 40 and 60% hydrates, water is not only adsorbed into the polymer matrix but also accumulated around the rolled chains.Such "free" water would not interact directly with the polymer chains but resemble bulk water.This conclusion agrees with the experimental results evidencing bulk water in the hydrates containing at least 30% of water by mass [52,55,56].

Dynamics of the Polymer Chains
The mean square displacement (MSD; cf.Equation ( 7)) is a measure of the polymer chains' mobility.Thus, the MSDs were calculated from the MD-generated trajectories.Since the polymers are not typical fluids, the subdiffusion occurs rather than Brownian motions.Consequently, the exponent α in the power law:

Dynamics of the Polymer Chains
The mean square displacement (MSD; cf.Equation ( 7)) is a measure of the polymer chains' mobility.Thus, the MSDs were calculated from the MD-generated trajectories.Since the polymers are not typical fluids, the subdiffusion occurs rather than Brownian motions.Consequently, the exponent α in the power law: is smaller than 1 rather than α = 1 typical of the diffusion described by the Einstein and Smoluchowski equation.The exponents α for hydroxyl O H3 and N M atoms of R37 with different water content are plotted in Figure 16.They were calculated by the least squares fitting of the logarithmic form of Equation ( 1): where K α is the generalized diffusion coefficient.The time interval for the fit was ca.500 ns-long, starting from the thirtieth ps of the simulation.
Molecules 2023, 28, x FOR PEER REVIEW different water content are plotted in Figure 16.They were calculated by the least s fitting of the logarithmic form of Equation ( 1): where Kα is the generalized diffusion coefficient.The time interval for the fit was ns-long, starting from the thirtieth ps of the simulation.The polymer chains meander around equilibrium points in the medium of co structures where the movement of the molecules is impeded.Thus, the average jumping from one equilibrium point to another is infinite.The more mobile the mo are, the higher α is in Equations ( 1) and (2).For each water concentration, the NM a the MPC phosphorylcholine group is more mobile than the OH3 of the HEMA hy one.Like the previously analyzed RDF functions, this suggests high flexibility of th side chains (cf.Figures 3 and 4) that favors the formation of diffusion channels [57,5 side chains are more mobile in the hydrated polymers because water acts as a plas Small water molecules penetrate the space between chains that increases free v [59], cf. Figure 15.Consequently, the hydrated polymers would be less viscous.
Molecular models of Rouse [60] and Zimm [61] applied for polymers in sem and dilute solutions, respectively, provide insight into polymer dynamics.In the model, the polymer chain is represented by beads connected by springs.Eac interacts with its neighbors through elastic forces and undergoes friction forces fr surroundings.Zimm extended Rouse's approach by considering hydrod interactions that alter the beads' dynamics.A force applied to one chain segment d The polymer chains meander around equilibrium points in the medium of complex structures where the movement of the molecules is impeded.Thus, the average time of jumping from one equilibrium point to another is infinite.The more mobile the molecules are, the higher α is in Equations ( 1) and (2).For each water concentration, the N M atom of the MPC phosphorylcholine group is more mobile than the O H3 of the HEMA hydroxyl one.Like the previously analyzed RDF functions, this suggests high flexibility of the MPC side chains (cf.Figures 3 and 4) that favors the formation of diffusion channels [57,58].The side chains are more mobile in the hydrated polymers because water acts as a plasticizer.Small water molecules penetrate the space between chains that increases free volumes [59], cf. Figure 15.Consequently, the hydrated polymers would be less viscous.
Molecular models of Rouse [60] and Zimm [61] applied for polymers in semidilute and dilute solutions, respectively, provide insight into polymer dynamics.In the Rouse model, the polymer chain is represented by beads connected by springs.Each bead interacts with its neighbors through elastic forces and undergoes friction forces from the surroundings.Zimm extended Rouse's approach by considering hydrodynamic interactions that alter the beads' dynamics.A force applied to one chain segment distorts the velocity field in the rest of the fluid influencing the motion of other segments.The distortion slowly vanishes with increased distance.The MSD obeys the following scaling behavior [62]: where τ m is the relaxation time of a monomer, and τ is the longest relaxation time of the chain.In the Rouse model, α = 0.5 for ideal chains and α = 0.54 for self-avoiding ones, while α = 0.67 in the Zimm model.Thus, the hydrodynamic interactions result in higher values of α.
The HEMA chain dynamics can be explained in terms of the Rouse model.In this case, there is no need to consider the hydrodynamic interactions in dry and hydrated polymers.MPC chains behave similarly in systems containing up to 20% water by mass independently of the simulation time and in those with 40 and 60% water until ca.30 ps is reached.At longer simulation times, α ≥ 0.6 suggests hydrodynamic interactions considered in the Zimm model.Thus, the simulated polymers are subjected to reptational movements, and interactions between neighboring side chains hamper their motions.Only in highly hydrated polymers can the MPC side chains move more freely thanks to the surrounding water molecules.

Dynamics of Water Molecules
Exponents α in Equation ( 1) for the MSD of water molecules in hydrated polymers were calculated in the same manner as those for O H3 and N M atoms.The only difference was the fitting range, which covered much longer time intervals exceeding 1000 ps.The values of exponents α reported in Table 4 suggest subdiffusion of water in hydrated copolymers at least up to the water content of 20% by mass, cf.Equation (3).In the systems containing 40 and 60% water, normal diffusion probably occurs as the exponents α are equal to 1 within the uncertainty limits of the fit at the confidence level of 95%.Higher exponents for the MPC-rich copolymers evidence that water is more mobile in these systems.The block or random arrangement of mers is not essential for mobility.These results are consistent with those discussed earlier.The hydroxyl and carbonyl groups of HEMA, and the carbonyl, choline, and phosphate groups of MPC form hydrogen bonds with surrounding water molecules and restrict their motions.In this manner, the so-called bound water arises.In highly hydrated systems (40 and 60% of H 2 O), water may resemble the bulk one to some extent.This is evidenced by the diffusion coefficients D (see Equation ( 8)) reported in Table 5.Their values amount from 53 to 58% of 2.84 × 10 −5 cm 2 s −1 for bulk water at 307 K [63].Thus, the decrease in mobility of water molecules caused by the polymer network is similar to that due to the cooling of water from 307 K to 285 K.In the latter temperature, D H2O = 1.62 × 10 −5 cm 2 s −1 [63].Such a phenomenon was also reported by other researchers [37,64].Although the diffusion coefficients are higher in copolymers with 37 MPC mers than in those with 13, the difference is small.Thus, flexible MPC chains can create efficient diffusion channels even in systems where they are outnumbered by HEMA mers.Moreover, water molecules in the outer layers of thick hydration shells of phosphorylcholine groups are fairly mobile [21][22][23][24][25][26].In the suggested two-state approach, the molecules' mobility in the bound water differs from that in the bulk-like one.Analysis of the self-term of the van Hove spacetime correlation functions (Equation ( 11)) supports this idea.Figure 17 shows G s → r , t for oxygen atoms of water molecules (O W ) in hydrated R37 in time intervals t from 1 to 100 ps long.For short simulation times, the distribution is unimodal.The single peak at r = 1.5 Å is due to vibrations of O W within local "cages" formed by the neighboring atoms.Initially, all atoms vibrate in their original positions, and only after some time can jump to the adjacent cage.Therefore, multiple peaks, at least two, suggest a "hopping" of O W from one cage to another for the longer simulation times.This is evident for the systems with 10 and 20% water.For more hydrated systems, the effect of "hopping" is masked by that of the self-diffusion of free water.Consequently, the self-correlation function is rather flat with a single maximum at r much bigger than the hydrogen bond length.
, x FOR PEER REVIEW 17 of 23 correlation functions (Equation ( 11)) supports this idea.Figure 17 shows  ( ⃗, ) for oxygen atoms of water molecules (OW) in hydrated R37 in time intervals t from 1 to 100 ps long.For short simulation times, the distribution is unimodal.The single peak at r = 1.5 Å is due to vibrations of OW within local "cages" formed by the neighboring atoms.Initially, all atoms vibrate in their original positions, and only after some time can jump to the adjacent cage.Therefore, multiple peaks, at least two, suggest a "hopping" of OW from one cage to another for the longer simulation times.This is evident for the systems with 10 and 20% water.For more hydrated systems, the effect of "hopping" is masked by that of the self-diffusion of free water.Consequently, the self-correlation function is rather flat with a single maximum at r much bigger than the hydrogen bond length.

Model Construction and Simulation Details
Molecular dynamics (MD) simulations were performed with the Amorphous Cell and Forcite modules of the Materials Studio package [65].To avoid calculation errors

Model Construction and Simulation Details
Molecular dynamics (MD) simulations were performed with the Amorphous Cell and Forcite modules of the Materials Studio package [65].To avoid calculation errors resulting from incorrect modeling of the starting structure, three modeled systems were constructed for each polymer configuration, and then the obtained results were averaged for each atomic system In the first step, the geometries of four polymer chains were optimized with the Smart method implemented in the Forcite module.The DREIDING force field (FF) [66] suitable for polymer systems [50,51,67,68] was applied.The main reason for choosing the DREIDING FF was the intention to model hydrogen interactions and hydrogen bonds.The proposed FF describes the hydrogen interactions explicitly by Lennard-Jones 12-10 potential [66].In this FF, the total potential energy E total is the sum of the interaction energy terms: The bond interaction energy E bond is due to chemical bonding.It comprises terms resulting from the bond stretching and bending (E B ), flat (E A ), torsion (E T ), and improper (E I ) angle components.The non-bond term (E non-bond ) encompasses van der Waals (E vdW ) and electrostatic (E Q ) interactions, and hydrogen bonds (E hb ) [66].
The maximum number of geometry optimization cycles and the energy convergence criterion were set to 1000 and 0.001 kcal/mol, respectively.Each isothermal-isochoric (NVT) MD simulation was 5-ns-long with a time step of 1 fs.To check the sufficiency of the simulation time, the non-hydrated B13, B37, R13, and R37 systems were simulated for 10 ns.In all cases, the system was equilibrated thermodynamically after 3 ns.The total energy of the system was kept constant for simulations longer than 3 ns.All structures reached thermodynamic equilibrium during the simulation.The equations of motion were integrated using the Verlet algorithm [69].The MD simulations were controlled by the Nosé-Hoover thermostat [70] set to a temperature of 307 K equal to the average for the eyeball surface [71].Copolymer configurations obtained in this manner were initial for the MD simulations of hydrates.
In the second step, simulations dealt with the polymers hydrated up to 60% swelling ratio.The initial hydrated systems were built with the Amorphous Cell module.To avoid long calculations, a single polymer chain occupied each unit cell.The sizes of the simulation unit cells adjusted to the final density of the system equal to 1.1 g/cm 3 are reported in Table 1.Methods of calculation were similar to those applied in the first step, except that the electrostatic and van der Waals interactions were calculated by the 3D Ewald summation.

Characteristics of the Simulated System
Several functions and quantities were obtained from the simulation results.Short definitions and relationships are given below for readers' convenience.More on this subject can be found in textbooks and monographs, e.g., in [72].
The radial distribution function (RDF) for atom pairs is defined conventionally as where A and B denote the atoms, ρ(r) is the number density of B at a distance r from A, and ρ 0 is the average density of B in the system.Hydrogen bonds D-H. ..A in the systems were identified based on the geometric criteria: the DA distance up to 3.5 Å, the HA distance shorter than 2.5 Å, and the DHA angle deviation from the straight one up to 30 • [48].The latter criterion is of marginal importance [73].
Free volume, V f , is operationally defined as that of the polymer material which is not occupied by the material itself but remains available for other molecules [74].To estimate the free volumes, the Connolly algorithm was applied [75].In this approach, the accessible molecular surface (the Connolly surface) is formed by virtually "rolling" a spherical probe over the van der Waals surface of a molecule, as illustrated in Figure 18.The free volume is outside the Connolly surface.
Hydrogen bonds D-H…A in the systems were identified based on the geometric criteria: the DA distance up to 3.5 Å, the HA distance shorter than 2.5 Å, and the DHA angle deviation from the straight one up to 30° [48].The latter criterion is of marginal importance [73].
Free volume, Vf, is operationally defined as that of the polymer material which is not occupied by the material itself but remains available for other molecules [74].To estimate the free volumes, the Connolly algorithm was applied [75].In this approach, the accessible molecular surface (the Connolly surface) is formed by virtually "rolling" a spherical probe over the van der Waals surface of a molecule, as illustrated in Figure 18.The free volume is outside the Connolly surface.Fractional free volume (FFV), f v , is the ratio of the free volume V f to the total volume of the system V: The following van der Waals radii for the atoms of the polymer molecules were applied in the calculations: r C = 1.70 Å, r H = 1.20 Å, r O = 1.52 Å, r N = 1.55 Å, and r P = 1.80 Å.
Mean square displacement (MSD), r 2 (t) , is defined as: where → r i (0) and → r i (t) are the initial and final positions of the i-th particle over the time interval t, and N is the number of equivalent particles.
The diffusion coefficient, D, is related to the MSD by Einstein's relation: The van Hove space-time correlation function, G → r , t , is defined as the correlation between the particle position → r + → r at the time t + t and its position → r at the time t .For a homogeneous system at equilibrium and the initial time t = 0, the G → r , t function depends on the relative distance: where the angular brackets represent the time average over the system and δ is the Dirac The self-term G s → r , t gives the distribution of displacements at time t.In this study, the G s → r , t functions for water molecules in the hydrated polymers were analyzed.

Conclusions
All simulated copolymers, dry and hydrated, proved to be amorphous.The backbone chains rolled inward with outward-facing side chains.Neither the copolymer composition nor the mers arrangement affected the HEMA sidechains' geometry.The latter remained fairly rigid even when the polymer was hydrated.Contrary to HEMA sidechains, those of MPC were flexible and coiled.That resulted in the mobility of the phosphorylcholine groups.However, the neighboring HEMA could stiffen the MPC chain.In hydrated polymers, particularly the two MPC-rich ones, water molecules penetrated voids between the MPC sidechains and uncoiled them to some extent.
Only interatomic distances suggested hydrogen bonds, as the molecules were not polarizable in these simulations.The bonds occur between the terminal hydroxyl groups of the HEMA sidechains and water molecules.The main polymer chain keeps water molecules from reaching the carbonyl oxygen atom of HEMA.Consequently, it weakens its propensity to hydrogen bonding.Nevertheless, the tendency is greater than that of the ether oxygen atoms.In this respect, carbonyl and ether groups of MPC sidechains did not differ from those of HEMA.However, the MPC proved more hydrophilic than HEMA due to the phosphate group being capable of forming hydrogen bonds through negatively charged oxygen atoms.Water molecules gathered around this and trimethylammonium groups.For the latter, this is due to electrostatic attraction only.Hydrogen bonding between water and nitrogen atoms is impossible because of the shielding methyl groups.
Non-hydrated block copolymers showed smaller free volumes than random ones.The fraction of volume available for water increased in the following order of the polymers: B13 < B37 < R13 ≈ R37.Thus, if HEMA mers separated the long MPC chains from one another, then penetration of water within the interchain spaces would be easier.Swelling of the polymers that accompanied hydration diminished the differences in the free volumes completely.
Flexible MPC side chains favored the formation of diffusion channels.A similar effect was reported for other flexible-chain polymers [57,58].Subdiffusion of water in hydrated copolymers was evident at least up to the water content of 20% by mass.The process consisted of the "hopping" of water molecules from one "cage" in the polymer matrix to another.Between the "hops", the molecule oscillates in the "cage" which is determined by the minimum of the potential energy of interactions.Further hydration increased the share of water resembling the bulk phase.Consequently, normal diffusion predominates in hydrated polymers with 40 and 60% of water by mass.Experimental studies confirmed such normal diffusion in polymers containing at least 30% of water [76].Moreover, water acted as a plasticizer increasing the mobility of polymer chains.
A general agreement of the conclusions from the present simulations with experimental results reported in the literature confirmed that molecular dynamics could facilitate the selection of potential materials for SCLs.

Molecules 2023 , 23 Figure 1 .
Figure 1.A fragment of the P(MPC-co-HEMA) copolymer chain.The atoms are tagged to simplify referring.

Figure 1 .
Figure 1.A fragment of the P(MPC-co-HEMA) copolymer chain.The atoms are tagged to simplify referring.

Figure 4 .
Figure 4. Snapshots of the simulated adjacent MPC side chains in non-hydrated block copolymer (a) and that containing 60% of water by mass (b).Hydrogen atoms and water molecules were omitted for the picture clarity.Atoms: C-grey, O-red, N-blue, and P-magenta.

Figure 3 .
Figure 3. Radial distribution functions g(r) of the C M -N M pair of MPC in the B13, B37, R13, and R37 copolymers.Blue line-non-hydrated polymer; black, red, green, and magenta lines-hydrated polymers containing 10, 20, 40, and 60% of water by mass, respectively.

Figure 4 .
Figure 4. Snapshots of the simulated adjacent MPC side chains in non-hydrated block copolymer (a) and that containing 60% of water by mass (b).Hydrogen atoms and water molecules were omitted for the picture clarity.Atoms: C-grey, O-red, N-blue, and P-magenta.

Figure 4 .
Figure 4. Snapshots of the simulated adjacent MPC side chains in non-hydrated block copolymer (a) and that containing 60% of water by mass (b).Hydrogen atoms and water molecules were omitted for the picture clarity.Atoms: C-grey, O-red, N-blue, and P-magenta.

Figure 5 .
Figure 5. Radial distribution functions g(r) for atom pairs: water oxygen (OW) and CH, OH3, and OH2 of HEMA for the R37 copolymer with various water content.Black, red, green, and magenta lineshydrated polymers containing 10, 20, 40, and 60% of water by mass, respectively.

Figure 5 .
Figure 5. Radial distribution functions g(r) for atom pairs: water oxygen (O W ) and C H , O H3 , and O H2 of HEMA for the R37 copolymer with various water content.Black, red, green, and magenta lines-hydrated polymers containing 10, 20, 40, and 60% of water by mass, respectively.

Figure 6
shows the RDF functions for MPC + water: C M -O W , N M -O W, and P M -O W pairs.The shortest distance between P M and O W , 4.1 Å, suggests that water molecules gather in the vicinity of the negatively charged phosphate group.Similar gathering around the positive trimethylammonium group also occurs, but at a longer distance of 4.7 Å because of the three blocking methyl groups.The distance between C M and O W , 5.0 Å, is equal to that between C H and O W . Indeed, the C H and C M atoms belong to the main polymer chain, and the HEMA and MPC side chains are identical up to the second O atom (O H3 and O M3 ).

Figure 6 .
Figure 6.Radial distribution functions g(r) for atom pairs: water oxygen (Ow) and CM, NM, and PM of MPC for the R37 copolymer with various water content.Black, red, green, and magenta lineshydrated polymers containing 10, 20, 40, and 60% of water by mass, respectively.

Figure 6 .
Figure 6.Radial distribution functions g(r) for atom pairs: water oxygen (O w ) and C M , N M , and P M of MPC for the R37 copolymer with various water content.Black, red, green, and magenta lines-hydrated polymers containing 10, 20, 40, and 60% of water by mass, respectively.

Figure 7 .
Figure 7. Radial distribution functions g(r) for atom pairs: water hydrogen and HEMA oxygen for the R37 copolymer with various water content.Green lines: OH1-HW, red: OH2-HW, and black: OH3-HW.

Figure 7 . 23 Figure 8 .
Figure 7. Radial distribution functions g(r) for atom pairs: water hydrogen and HEMA oxygen for the R37 copolymer with various water content.Green lines: O H1 -H W , red: O H2 -H W , and black: O H3 -H W . Molecules 2023, 28, x FOR PEER REVIEW 9 of 23

Figure 8 .
Figure 8. Radial distribution functions g(r) for atom pairs: water oxygen and HEMA oxygen for the R37 copolymer with various water content.Green lines: O H1 -O W , red: O H2 -O W , and black: O H3 -O W .

Figure 9 .
Figure 9. Radial distribution functions g(r) for atom pairs: water hydrogen and MPC oxygen for the R37 copolymer with various water content.Black lines: O M1 -H W , red: O M2 -H W ; green: O M3 -H W , blue: O M(4,5) -H W , cyan: O M6 -H W . Answering the second question requires a more detailed discussion.The RDFs for O M4 and O M5 are pretty much the same because the two atoms are equivalent in an MPC molecule, as illustrated by the resonance structures:

Figure 10 .
Figure 10.Radial distribution functions g(r) for atom pairs: water hydrogen and MPC nitrogen or phosphorus for the R37 copolymer with various water content.Black lines: NM-HW, red: PM-HW.

Figure 10 . 23 Figure 11 .
Figure 10.Radial distribution functions g(r) for atom pairs: water hydrogen and MPC nitrogen or phosphorus for the R37 copolymer with various water content.Black lines: N M -H W , red: P M -H W . Molecules 2023, 28, x FOR PEER REVIEW 11 of 23

Figure 11 .
Figure 11.Radial distribution functions g(r) for atom pairs: water oxygen and MPC oxygen for the R37 copolymer with various water content.Black lines: O M1 -O W , red: O M2 -O W ; green: O M3 -O W , blue: O M(4,5) -O W , cyan: O M6 -O W .

Figure 12 .
Figure 12.Radial distribution functions g(r) for atom pairs: water oxygen and MPC nitrogen or phosphorus for the R37 copolymer with various water content.Black lines: NM-OW, red: PM-OW.

Figure 12 .
Figure 12.Radial distribution functions g(r) for atom pairs: water oxygen and MPC nitrogen or phosphorus for the R37 copolymer with various water content.Black lines: N M -O W , red: P M -O W .

Figure 13 .
Figure 13.Morphology of free volumes in R37, dry (a) and with different water content: 10% (b), 20% (c), 40% (d), and 60% (e) by mass for the probe of infinitesimal radius.The walls of the empty channels making up the free volumes are marked in blue.

Figure 13 .
Figure 13.Morphology of free volumes in R37, dry (a) and with different water content: 10% (b), 20% (c), 40% (d), and 60% (e) by mass for the probe of infinitesimal radius.The walls of the empty channels making up the free volumes are marked in blue.

Figure 14 .
Figure 14.The fractional free volume (FFV) of non-hydrated P(MPC-co-HEMA) as determined with probes of various radii.

Figure 14 .
Figure 14.The fractional free volume (FFV) of non-hydrated P(MPC-co-HEMA) as determined with probes of various radii.

Figure 15 .
Figure 15.The volume in the simulated R37 material available for or filled with water: dry (a) and hydrated copolymers containing 10% (b), 20% (c), 40% (d), and 60% (e) of water by mass.The walls of polymer channels available for the probe with a radius of 1.4 Å (i.e., a "water molecule") are marked in blue.

Figure 15 .
Figure 15.The volume in the simulated R37 material available for or filled with water: dry (a) and hydrated copolymers containing 10% (b), 20% (c), 40% (d), and 60% (e) of water by mass.The walls of polymer channels available for the probe with a radius of 1.4 Å (i.e., a "water molecule") are marked in blue.

Figure 16 .
Figure 16.The exponents in the power law (Equation (1)) that evidence subdiffusion of (open symbols) and NM atoms (filled symbols) of R37 vs. the mass concentration of water.L guides for the eye only.

Figure 16 .
Figure 16.The exponents in the power law (Equation (1)) that evidence subdiffusion of the O H3 (open symbols) and N M atoms (filled symbols) of R37 vs. the mass concentration of water.Lines are guides for the eye only.

Figure 17 .
Figure 17.Van Hove self-correlation functions G s → r , t calculated for water molecules in R37.

Figure 18 .
Figure 18."Rolling" of a spherical probe with radius r over the surface of a molecule built of two atoms with van der Waals radii r 1 and r 2 .
delta function.The G → r , t function can be split into the self-term G s → r , t for i = j and distinct term G d → r , t for i = j: G → r , t = G s → r , t + G d → r , t ,

Table 1 .
The cell sizes for the P(MPC-co-HEMA) systems after 5-ns-long MD simulations.

Table 1 .
The cell sizes for the P(MPC-co-HEMA) systems after 5-ns-long MD simulations.

Table 2 .
The fractional free volumes (FFVs) in P(MPC-co-HEMA) of different built and degrees of hydration for the probe of infinitesimal radius.

Table 3 .
The fractional free volume (FFVs) in P(MPC-co-HEMA) available for water molecules, calculated with a spherical probe of a radius of 1.4 Å.

Table 3 .
The fractional free volume (FFVs) in P(MPC-co-HEMA) available for water molecules, calculated with a spherical probe of a radius of 1.4 Å.

Table 5 .
The diffusion coefficient of water in hydrated P(MPC-co-HEMA) copolymers.