Modelling thermoresponsive polymer brush by mesoscale computer simulations

We consider a functional surface comprising thermoresponsive polymer chains, the material that has found numerous technological and biomedical applications. However, to achieve the required time and length scales for computer modelling of such applications, one is compelled to use coarse-grained mesoscopic modelling approaches. The model used here is based on the previous work [Soto-Figueroa et al., Soft Matter, 8, 1871 (2012)], and it mimics the principal feature of the poly(N-iso-propylacrylamide) (PNIPAM), namely, the rapid change of its hydrophilicity at the lower critical solution temperature (LCST). For the case of an isolated chain, we discuss scaling properties of the radius of gyration, end-to-end distance, various distribution functions, and the density profile of monomers below and above the LCST. For the case of the model thermoresposive brush, we search for the optimum grafting density at which the change in the brush height, upon crossing the LCST, reaches its maximum value. The interpretation of the thermoresponse, in terms of the Alexander-de Gennes blobs and the level of solvation of polymer chains in a brush, is provided.


Introduction
Stimuli-responsive polymers have a wide range of applications in biotechnology and medicine [1]. Thermoresponsive polymers form a sub-class of these materials, and the PNIPAM is their prominent representative. Applications of this polymer are based on a drastic change of its solubility in water that occurs around the LCST. Below LCST, PNIPAM is hydrophilic and it swells in water, whereas above LCST, PNIPAM chains collapse [2]. The LCST for this polymer is close to 32 • C, which is important from the point of view of the physiological applications [3].
Stability of a swollen state of PNIPAM below LCST is attributed to the formation of water cages and hydrogen bond bridges with water molecules [4,5]. Upon the temperature increase, the number of PNIPAM-water hydrogen bonds decrease, whereas the number of the polymer-polymer hydrogen bonds increase, accompanied by desolvation of the hydrophobic groups, as suggested by spectroscopy studies [6][7][8][9]. The desolvation of the polymer above the LCST manifests itself in an essential decrease in the hydration number [10,11].
This microscopic picture is confirmed in a number of computer simulation studies of a single PNIPAM chain using a fully atomistic modelling approach [12][13][14]. It was shown that a network of hydrogen bonds is formed between the PNIPAM segments and the water molecules below LCST resulting in hydration of the PNIPAM chains and their coiled conformations. Above LCST, such a network breaks in the proximity of hydrophobic isopropyl groups, and the PNIPAM chains collapse into compact globules [12]. The transition, that occurs at LCST, can therefore be classified as the conformational coil-to-globule transition. Dalgicdir et al. [13] tested the OPLS/AA force field for a single PNIPAM 40-mer in water. They found that the experimental coil-to-globule collapse enthalpy of PNIPAM in water can be reproduced by simulations forcing the amide-water electrostatic interactions stronger than in the original OPLS model. in the form of a patterned mixed brush, with one of the components being made of the PNIPAM chains of considerable length. Below LCST, the PNIPAM areas of a brush swell and rise up, breaking the proteins out free. Considering that the swelling of a small PNIPAM region of such brush can be described by the mid-level coarse-grained models discussed in a previous paragraph, the role of the spatial pattern in the effectiveness of such smart surface clearly requires large system sizes inaccessible to such models. Here, such higher level of coarse-graining, as introduced above, would be of a great aid.
Secondly, there is already some work done on modelling the polymers that involve the PNIPAM fragments [28][29][30][31]. Justification of the effective interaction parameters between the PNIPAM monomers and water, including the formation of hydrogen bonds, was discussed by Vicente with colleagues [29]. They addressed both theoretical calculation of Flory-Huggins parameter and the Monte Carlo simulations approach involving the evaluation of the coordination numbers and binding energies between interacting particles. As a result, the temperature dependence of the parameter for the PNIPAM chain was obtained.
The aim of this study is to examine whether or not this type of modelling faithfully reproduces the main features of the PNIPAM below as well as above the LCST. Both the cases of an isolated PNIPAM chain and the ensemble of the PNIPAM chains, arranged as a polymer brush, are considered. For the former, we examine the scaling behaviour of a number of its properties with respect to the chain length. For the latter, we focus on a search for an optimum grafting density of a brush at which the thermoswitching effect has the maximum magnitude, and also discuss microscopic mechanisms behind the thermoswitching behaviour.
This research work is dedicated to Taras Bryk, leading expert in the dynamic properties of liquids, liquid metals and alloys, a specialist in the ab initio and classical molecular dynamics simulations, Director of the Institute for Condensed Matter Physics, on the occasion of his 60th birthday. We wish him scientific prolificity, new ideas and achievements for the years to come.
The outline of the study is as follows. In section 2 we provide details of the model and of the simulation technique, the properties of an isolated PNIPAM chain are discussed in section 3, whereas section 4 contains the results for the case of the PNIPAM brush at different grafting densities. The study is summarised by conclusions.

Modelling details
We use the dissipative particle dynamics (DPD) simulation approach, which belongs to the class of mesoscopic simulation techniques. In doing this we closely follow the seminal work by Groot and Warren [32]. In this approach, the monomers of a polymer chain, as well as the solvent, are all represented as soft spheres of equal size that interact via pairwise forces. The diameter of the soft bead defines a native length-scale in a model, and the energy unit is chosen equal to B = 1, where B is the Boltzmann constant and is the temperature, time unit is * = 1. The monomers are bonded via harmonic springs which results in the bonding force: where x = x − x , x and x are the positions of the -th and -th bead, and is the associated energy. The total non-bonded force F acting on -th bead from its -th counterpart can be expressed as a sum of three contributions: where F C is the conservative force caused by the repulsion of the -th and -th beads, F D is the dissipative force responsible for friction between beads (similar to the Stokes force), and the random force F R which, together with the dissipative force, maintains a constant temperature of the system. Expressions for these terms are pretty standard and are given below: 33302-3 , v i is the velocity of the -th monomer, and is the amplitude for the conservative repulsive force. The dissipative force has an amplitude and its decrease with the particle separation is given by the function D ( ). Similarly, the amplitude for the random force is and its decrease is provided by R ( ).
is a Gaussian random variable. As shown in [33], to satisfy the requirement of a detailed balance, the relation 2 = 2 and D ( ) = [ R ( )] 2 should hold. Here, we use a quadratic decaying form for the weight function: Mesoscopic modelling of the PNIPAM polymer follows closely the work by Soto-Figueroa et al. [28], where the PNIPAM-containing block-copolymers were considered. Namely, each repeating unit of the PNIPAM is represented via a single soft spherical bead, as shown in figure 1. Each polymer chain consists of such monomers that are bonded pairwise via the force (1) with the spring constant = 4. In general, repulsive strength of the conservative force, , in equation (3) depends on the temperature and density. In their practical implementation of the DPD approach, Groot and Warren considered the fixed bulk number density of beads equal to 3, and parameter was found from matching the inverse compressibility of a model to that of water at normal conditions [32]. This resulted in the value of = 25. We denote polymer beads by a subscript " ", and those representing water, by " ". It is assumed that for the polymer-polymer and water-water interaction, one has = = 25. Mixed, polymer-water interaction strength, , depends on the polymer level of hydrophilicity. In the case of PNIPAM, one expects ≈ 25 below LCST and > 25 above it. Hence, the strength of the polymer-water hydrogen bonding does enter this model implicitly, via the value of the . The temperature dependence of was estimated from that for the Flory-Huggins mixing parameter , see, figure 3 in [28], using the relation between the two [32]. To simplify the analysis, in this study we concentrate on two temperatures only: one below, and another one above the LCST. Namely, at = 298 K (below LCST), the estimate for is ≈ 25.6, whereas at = 310 K (above LCST), one has ≈ 38. The simulation setup is as follows. The simulation box dimensions are , and , which are specified in each case being considered. The periodic boundary conditions are applied along both OX and OY axes, whereas two walls representing a substrate are imposed at = 0 and = . Monomer-surface interaction plays a big role in specifying the structure of a brush at various densities [34]. Here, we restricted our attention to a singular case of repulsive impenetrable walls. The repulsion from the wall is of the same soft nature as between the beads and is realized using the expression (3). Here, x denotes the vector that runs from the center of -th bead towards the nearest wall and is normal to it. The amplitude of the bead-wall repulsion is fixed at = 25. This type of repulsion does not ensure the impenetrability of the walls. Therefore, the elastic reflection is imposed when the bead crosses the either of the walls.
In general, we consider a certain number of chains, and one of the end beads of each chain is grafted to a specified fixed point of the wall located at the = 0 plane. Grafting is performed by using 33302-4 the harmonic force (1) with = 4. If grafting points are chosen randomly, then one ought to average over many such arrangements and, additionally, the results are affected by local inhomogeneities of the surface density of grafting points. To reduce these effects, we opted to arrange the grafting points on a regular grid and perform a single simulation run for each . Specifically, the grid of × grafting points is made in the form of a square lattice, where is the integer value of √ . The rest − 2 grafting points are distributed randomly across the wall. The grafting density is measured as the surface number density of grafts, = /( ). The chains are not allowed to slide on the surface. The rest of the simulation box is filled-in with solvent beads to achieve a required total bulk density of beads equal to 3 [32]. The time step was chosen equal to Δ = 0.04 in the DPD units. For the case of an isolated chain, = 1, considered in section 3, 8 · 10 6 DPD steps are performed, where the first 10 6 of them are discarded to allow conformational relaxation of the chain. For the case of a brush, > 1, covered in section 4, 10 6 DPD steps were used, where first 3 · 10 5 of them were discarded for relaxation purpose.

Isolated grafted thermoresponsive chain
Scaling properties of polymer chains, initiated in the seminal works by Flory, des Cloizeaux, de Gennes [35][36][37], are considered to be a well studied topic now. These studies brought the concept of universality into the field of polymer physics, mirroring a self-similarity of polymeric conformations. The universality features hold especially well at a mesoscopic level, where the fine chemical details are smeared out because of coarse-graining. Some examples can be found in our previous studies [38][39][40][41][42]. As was discussed in section 2, the PNIPAM polymer can be mapped onto a linear chain with the polymerwater interactions specified by the magnitude of in the expression of a soft repulsive force (3). Below LCST, the value of ≈ 25.6 is found [28], which is pretty close to the case of athermal solvent, = = = 25. Such case was extensively studied before, see, e.g. [38,39] and references therein. The aim of this study is not to merely repeat these findings but to perform a detailed comparison between the properties of an isolated chain at = 298 K (below LCST) and at 310 K (above it). Such comparison, that is seen as the first step in studying the thermoresponsive brush to be performed in section 4, was not done in detail before.
We consider an isolated chain that is grafted to the plane = 0 and surrounded by a solvent that represents water. The details of modelling are provided in section 2. The set of chain lengths, = 10, 14, 20, 28, 40, and 56, are considered in the three dimensional space. The simulation box dimensions in each case is chosen: (i) to be large enough to accommodate the most probable conformations of a chain, and (ii) not to be oversized to save on a computer time. Some recipes were suggested earlier [39], and in the current study we found a good compromise by the following choice: [38] is an equilibrium bond length between the monomers for the chosen values of and in equations (1) and (3) and a bulk number density of 3. In parallel, we perform similar simulations for a free chain of the same length at the same solvent conditions. The periodic boundary conditions are applied along all three spatial axes in this case. Simulation runs are of the duration of 8 · 10 6 DPD steps for each chain length , and the first 10 6 steps are discarded from the analysis to allow equilibration of the system. Let us consider the size properties of a polymer chain. Most generally, the spatial distribution of its mass is given by the gyration tensor, G, with its components defined as [43] Here, , and denote Cartesian coordinates of -th monomer and these for the centre of mass of a chain, respectively, and = 1, 2, 3.
The eigenvalues of the gyration tensor are sorted as follows 1 > 2 > 3 . These allow one to evaluate the instantaneous squared radius of gyration of a chain, 2 = 3

=1
. The average value for the radius of gyration is defined as follows:  where the averaging is performed over the time trajectory and we already used the scaling law with respect to the number of monomers, , [36,37]. Here, is a pre-factor, and is the scaling exponent which depends on the solvent quality. In particular, ≈ 0.5882 [44] for the case of good solvent, = 1/2 for the random walk (Gaussian chain) and = 1/3 for the case of poor solvent. Grafted chain is expected to obey the power law with the same exponent as its free counterpart but with a different pre-factor [45].
Simulation data obtained for¯and their numeric fits to the from (8) are shown in figure 2. In each case, the data are displayed in a log-log scale to allow clear visualisation of the quality of the fits represented as straight lines. The first obvious thing is that¯is essentially higher at = 298 K than at 310 K (where the chain is in a state of collapse). This will have important consequences for the case of a polymer brush, as discussed later, in section 4. All data shown in figure 2 are found to fit rather well to the scaling law (8), providing a robust estimate for the critical exponent . The estimates for ≈ 0.582 − 0.586 at = 298 K are in a very good agreement with the known most accurate result ≈ 0.5882 [44] for the Flory exponent. At = 310 K, however, the value of ≈ 0.300 is close to but lower than the expected exact result, = 1/3, indicating a strong mutual penetration of monomers because of their soft nature. Grafting a chain to a surface indeed affects the value of a pre-factor but the difference is rather minor, i.e., within the accuracy of simulations.  Similar results are obtained for the average end-to-end distance of a chain defined viā pre-factor and the critical exponent have the same meaning as in equation (8). In general, the end-toend distance is found to be more demanding in terms of the simulation length and it displays a scaling behaviour starting from longer chain lengths, ⩾ 20. The estimates for are acceptable for the case of = 298 K, whereas the value ≈ 0.27 at = 310 K is even lower than that obtained from the analysis of¯.
as shown in figure 4. The function ( ) is suggested by Lhuillier [46], following his considerations on the probabilities of fully collapsed and fully stretched conformations. Here, ≈ 42 is the pre-factor, the first term within the parenthesis is dominant at small * , whereas the second term defines the rate at which ( * ) decays at large * . We also would like to point out some details of our previous analysis of this distribution for the case of a free chain [38,39].
Scaled distributions are also examined for the case of a poor solvent, = 310 K. We observed that the initial distributions, in terms of , shift towards larger and increase moderately in their height upon the increase of the chain length, . Therefore, the following expression for the shifted radius of gyration, ′ , is proposed in this case (note that in the DPD simulations the length unit is assumed to be equal to unity)  figure 4 and one may note a rather good match of the simulation data for different . The respective master curves at both temperatures, = 298 K and 310 K, are remarkably similar in their shape, but we did not explore this similarity in any deeper terms here. Now, we discuss similar distributions for the end-to-end distance, . Here, only the case of a good solvent, = 298 K, is considered. The dimensionless end-to-end distance, * = /¯, is introduced, where the scaling laws for¯are shown in the left-hand frame of figure 3. Similarly to the case of the  , scaling of the distributions ( * ), obtained at different chain length , enables us to introduce the master curve see figure 5. The analytic form for the function ( ) is constructed following the analysis by des Cloizeaux [47] and de Gennes [48]. Using numeric fits for the simulation data, we found that ≈ 5.65, ≈ 2.63 and ≈ 1.15. Since > 1, the shape of ( * ) at small * is convex. One should mention that the theoretical analysis [47,48] yields the value of = (1 − )/ ≈ 0.27 < 1 resulting in a concave curve at small * . We believe that this discrepancy is the result of the excluded volume effects at small in the DPD modelling approach, and these reduce the probability for the values * → 0.  (13) and (14).
Finally, we examine the density profile ( ) for a single model PNIPAM chain grafted to the plane = 0, and compare it to its counterpart for a free chain. In the latter case, the reference plane is oriented randomly in space. It intersects the first bead of a chain and is completely transparent. For the case of = 298 K, the density profile can be built in terms of the dimensionless separation from the wall, * = /¯using the same scaling law for¯as above. The data for a scaled profile, ( * ), obtained at different , is found to very well follow their respective master curves (for the grafted and free chain, respectively), see the left-hand frame of figure 6.
In the case of a grafted chain, the plane = 0 is impenetrable and it causes depletion of the polymer beads in its vicinity, see a respective plot in the left-hand frame of figure 6. The analytic form for the scaled profile ( * ) can be constructed following considerations by de Gennes [37]. In particular, at small , one obtains ( * ) ∼ ( * ) 1− . At large * , we assume an exponential decay of ( * ) similarly

33302-8
to equation (10), but with unknown parameters. Then, we assume that both terms can be multiplied Numeric fits of simulation data yield ≈ 0.3, ≈ 0.56, and ≈ 1.74 and the obtained master curve represents simulation data very well, see respective plot in the left-hand frame of figure 6. In the case of a free chain, there is no effect from the impenetrable wall. Therefore, a profile has the Gauss-like shape, see a respective plot in the left-hand frame of figure 6. This leads to a guess that the analytic form for the scaled profile ( * ) in this case can be constructed by modifying the equation (13) by: (i) removing the term * 1− , and, (ii) shifting the argument in the exponential term Note, that the same coefficients, , , and the exponent enter this expression, whereas the amount of shift is given by ≈ 0.247. This analytic form very well describes the scaled profile ( * ) for a free chain, too, see the left-hand frame of figure 6.
The situation is quite similar at = 310 K, in which case the scaling law for¯, shown in the righthand frame of figure 2, is used. The DPD simulation data for both grafted and free chain, represented in terms of * = /¯, reasonably well follow their respective universal shapes ( * ), see right-hand frame of figure 6. One should remark that both shapes are very similar to their counterparts at = 298 K, see left-hand frame of the same figure, but their real separation from the wall is different because of a different scaling exponent for¯. Similarly to the case of figure 4, we have not performed analytic fits for the ( * ) distributions at = 310 K.
To conclude this section, here we performed validation of the scaling behaviour of a single chain representing the PNIPAM polymer at a mesoscopic level. Following reference [28], the temperature effect is modelled via the change of the repulsive strength between PNIPAM monomers and the beads representing water, see equation 3. Both cases of = 298 K (below the LCST), and = 310 K (above it) are considered. At the first stage of this analysis, the scaling laws for the average radius of gyration,¯, and these for the end-to-end distance,¯, are obtained by performing simulations at a range of chain lengths, = 10-56. The results are shown in figures 2 and 3 and these indicate a reasonably good agreement with available theoretical and computational results. At the following stage, we examined the scaling properties of the distributions obtained for the and values and the density profile of monomers normally to the wall. Where known, the asymptotic scaling behaviour was used from the available literature, otherwise, suitable analytic forms were constructed leading to a very good mapping of the obtained simulation data onto master curves. Therefore, the mesoscopic model for the PNIPAM polymer employed in this study, rather accurately demonstrates the scaling properties for the size properties,¯and¯, their distributions, and the density profile of monomers in relation to the wall.

Thermoresponsive polymer brush
We proceed now to the case of a polymer brush formed of grafted chains according to the procedure explained in detail in section 2. The length of all chains is fixed at = 60. The new parameter that is brought into a play here is the grafting number density of chains, , introduced in section 2. At low , polymer chains behave independently of one another exhibiting the so-called "mushroom" regime. Here, the behaviour of chains is principally governed by the quality of the solvent, the situation being similar to the case of an isolated chain discussed in section 3. At a very high , the chains are closely packed, and their properties are affected mostly by the effects of the excluded volume from the overlaps with the neighbouring chain. The crossover between these two regimes is expected at the intermediate values of . Scaling properties of polymer chains are widely covered in a number of studies [34,37,45,[49][50][51][52][53]. Therefore, here we restrict our analysis to the aspects most relevant to the thermoswitching behaviour of the model chains representing the PNIPAM polymer brush.
Let us introduce some relevant properties of a brush. The effective brush height can be defined as [52]

33302-9
where ( ) is the density profile for the brush beads, and is the simulation box dimension along the OZ axis. Besides this collective feature of the ensemble of chains, we also consider individual chain characteristics related to their conformations. The size characteristics such as¯and¯, are already introduced in equations (8) and (9). We also complement our analysis by the shape characteristics such as: the average asphericity, , and the shape descriptor, , [39,41,42] defined as , and the averaging in (16) is performed over the time trajectory and over all chains in the brush. The benefit from the use of the shape descriptor is that it is symmetric, namely = −1 for an infinitely thin disc, = 1 for infinitely thin rod and = 0 for the case of a sphere.
It is well known that with an increase of the grafting density a brush enters the so-called "dense brush regime", where individual chains are stretched out and the brush height increases [45]. To quantify the difference when and to what extent this takes place at = 298 K and 310 K, we consider the set of ratios † =¯( = 298 K) between respective values obtained at these two temperatures. The dependencies of these ratios on the grafting density are shown in figure 7. The effectiveness of the thermoresponsive brush as a functional surface can be measured via the magnitude of the ℎ † ratio, which is found to be the largest at ⩽ 0.3. One can also see that all the other fractions, † , † , † , and † , reach their maximum values at ≈ 0.3. This hints at this interval of grafting densities as the one where the difference between two temperatures in terms of the amount of a stretch of individual chains is the largest.
The narrow interval of grafting densities around ≈ 0.3 is a good candidate for designing a thermoswitchable brush from an additional reason as well. Figure 8 shows the structure of a brush surface at = 310 K, in the collapsed state. It can be seen that the surface is very inhomogeneous showing islands of collapsed chains and the regions of exposed surface until the grafting density reaches the value of ≈ 0.3. Such a behaviour at low grafting densities is well known since the early simulation studies [51,54]. In most practical applications, such inhomogeneity would be rather unwelcome as it may hamper the adsorption of nanoparticles in the collapsed state of a brush. Additionally, application of the PNIPAM-based brush as a part of a mixed brush system for thermocontrolled adsorption/desorption of proteins [27] requires the highest possible grafting density of the PNIPAM component to ensure strong push-effect at the desorption stage. These considerations lead us to believe that the grafting densities close to the ≈ 0.3 are optimum for the application of a thermoresponsive brush as an effective height-changing functional surface.
Let us examine the origins of the differences in the brush behaviour at various grafting densities. To this end, the blob interpretation of polymer chains by Alexander and de Gennes [37,55] can be invoked. In brief, if the average separation between the grafted points, , is getting smaller than the average radius of gyration of an isolated chain,¯, then the chain monomers started to redistribute their mass normally to the surface because of the excluded volume effect. The chain can be interpreted as a linear stack of blobs of diameter , assembled normally to the surface, where the number of monomers in one blob can be estimated from the relation = , where = 0.5882, 0.5, and 1/3 for the good solvent, -point and poor solvent regime, respectively, inside each blob. The blob interpretation leads to the following scaling law forh [37]h

33302-10
where ≈ 0.35 for the case of a good solvent within a blob, and = 1/2 for the case of a -point within a blob (desolvated chain). Similar scaling laws are also expected for both the values of¯and¯for each chain in a brush. Hence, the changes in the brush height upon increase of are governed by two factors. The first one is the excluded volume effect between the neighbouring chains, that is brought into play when ( ) decreases below the value. The dependence of ( ) on is easily established, ( ) = −1/2 , and it is shown in the left-hand frame of figure 9. As one can see, at = 298 K the value of drops below at ≈ 0.07, whereas at = 310 K this occurs for a much higher grafting density, ≈ 0.5. As the result, at grafting density = 0.3, the brush at = 298 K is already deep into the dense brush regime, whereas at = 310 K it is still in the "mushroom" state, thus explaining a large difference in its height  at these respective temperatures.
The second factor is the solvation regime of chain monomers within each blob, which affects the value of the exponent . To clarify this factor, we performed analysis of the contents of the first coordination sphere for all monomers at each grafting density. The radius of the first coordination sphere was chosen equal to the cutoff distance of the conservative force (3), and the number of polymeric monomers and that for solvent beads were counted. These were subsequently averaged over all existing monomers yielding and , the partial coordination numbers associated with the polymer-polymer and polymer-solvent pairs. Their dependencies on the grafting density are shown in the right-hand frame of figure 9. Total number of beads in the first coordination sphere, + , is close to six, as expected. At a low grafting density, one has strong solvation of polymer beads, ≫ , whereas strong desolvation occurs at high grafting densities, where ≪ . Crossover takes place at the grafting densities close to ≈ 0.45 (see, black dashed line in figure 9). Let us note that at = 0.3 (i.e., the supposed optimum grafting density for the thermoswitchable brush) the polymer chains are still, at least partially, solvated. Now, after the solvation regimes at various grafting densities are clarified, the scaling laws of the form (18) can be verified. To do so, we take into account that = 60 = const in current simulations, and assume the presence of some constant terms, ℎ 0 , 0 , and 0 ℎ , in the scaling expressions Then, we attempted to find the respective intervals for , where the fits with = 0.35 (solvation regime) and = 0.5 (desolvation regime) are reasonably accurate. The results are shown in figure 10. The following values for the coefficients entering equation (19) are found: ℎ 0 ≈ 2.9, ℎ ≈ 19.0, 0 ≈ 2.2, ≈ 3.6, 0 ≈ 4.2, and ≈ 13.4. The intervals being found are: 0.1 < < 0.5 for the = 0.35 fit and 0.5 < < 1.3 for the = 0.5 fit. These are in a good agreement with the picture of solvation regimes, as identified in figure 9. Therefore, we may conclude that, in general, the model polymer brush discussed here follow the scaling laws for the brush height, radius of gyration and the end-to-end distance, within the blob interpretation suggested by Alexander and de Gennes [34,37,45,55].
Therefore, computer simulations of the model thermoresponsive brush considered in this study reveals that there is some optimum grafting density, ≈ 0.3, at which the ratio between the brush height at = 298 K and 310 K is close to its maximum value (see, figure 7) and the structure of the brush in its collapsed state is homogeneous (see, figure 8). There is a combination of two factors that contribute to such a difference between the state of the brush for = 0.3 and at these two temperatures. The first  one is related to the excluded volume effect between the neighbouring chains for = 0.3. It already promotes the increase of the brush height at = 298 K [ ( ) <¯], but has no influence on it at 310 K [ ( ) >¯], see left-hand frame of figure 9. The second factor is related to the fact that the monomers are, at least partially, solvated for = 0.3 at = 298 K, see right frame of figure 9. This causes the chains to swell more as compared to the desolvated chains at = 310 K. With the increase of , the chains at = 298 K became desolvated as well, thus decreasing the difference between the cases of = 298 K and 310 K. With an increase of the grafting density towards dense brush regime, brush chains became equally desolvated and stretched both at = 298 K (below LCST) and = 310 K (above LCST), and the difference between these two temperatures vanishes, see, figure 7.

Conclusions
Thermoresponsive polymer brushes, based on the PNIPAM polymer, are widely used in a number of practical applications. These include drug delivery systems, protein and microalgae adsorption/desorption, membrane based filtering, etc. Their funtionality is based on the change of hydrophilicity of PNIPAM chains in water solution when the temperature crosses the LCST point. In particular, below LCST, PNIPAM is hydrophilic due to its capability to form hydrogen bonds with water molecules and, as a result, its chains adopt swollen coil-like conformation. Above the LCST, PNIPAM turns into a hydrophobic polymer and its chains collapse into a globule. Clear understanding of all the details of such coil-to-globule transition, covering the role of the chain length, details of molecular architecture, brush density and other factors, are of much benefit for designing new thermoresponsive functional materials.
Computer simulations are capable of shedding much light towards clarifying these aspects, but modelling the functionality of the thermoresponsive brush faces a number of challenges related to the need to cover relatively large time-and length-scales. In this study we employed the mesoscale simulation approach of the DPD with the use of the parametrisation for the PNIPAM chain that was performed earlier [28]. In this approach, the temperature dependent hydrophilicity is mapped onto the effective parameters for the polymer-solvent interaction. The principal aims of our study were to test the simple linear model for the PNIPAM polymer both below and above the LCST, and to get some insight towards the optimal brush structure to ensure its functioning as a thermoswitchable functional surface.
For the case of isolated grafted chain, we examined its scaling properties above and below the LCST and found the scaling laws to very well match the known results. In particular, the model reproduces the correct scaling behaviour for the radius of gyration, end-to-end distance, as well as distributions of their values and the density profile with respect to the distance to the wall.
For the case of a brush, the model reveals the presence of an optimum grafting density, ≈ 0.3, at which the ratio between the brush height below LCST and above it reaches its maximum value, indicating the largest amplitude of a thermoswitching effect. We note that this ratio is about 2.5, close to some reported experimental values [27]. At this grafting density, similar ratios, evaluated for the radius of gyration, end-to-end distance, asphericity, and the shape factor of individual chains, all reach their respective maxima. The reason for such a behaviour is examined by using the concept of blobs for a polymer brush by Alexander and de Gennes [37,55]. At a low grafting density, the brush is in a mushroom regime, i.e., it comprises a set of non-interacting chains. At a higher grafting density, the excluded volume effect between the neighbouring chains became important and causes the stretching of the polymer chain away from a wall leading to the interpretation of each chain as a linear stack of blobs. We found that at the grafting densities of ≈ 0.3, there is an essential difference between the chain conformations below and above the LCST. In particular, below LCST, the chains of a brush are deep in a stretched regime, whereas above LCST, the chains are still in a mushroom regime. This difference in chain conformations is further amplified by a good solvent condition below LCST and by poor solvent conditions above it. Both factors combined are responsible for the largest amplitude of a thermoswitching effect observed at this grafting density.
The results of this study can be used in a number of ways. First, they open up the possibility to model of the temperature controlled adsorption/desorption of proteins, as, e.g., in the experimental work reported in reference [27], or other macromolecules. The results can also be used to predict an optimal grafting density of a thermoresponsive polymer brush. For instance, one may perform the measurements 33302-13 of the radius of gyration of isolated chains first, and then use these to estimate the grafting density when the mushroom regime switches into the stretched one. Other applications of the presented results are possible as well.