# Human-Robot Interaction using VAHR:
Virtual Assistant, Human, and Robots in the Loop

Ahmad Amine1, Mostafa Aldilati1, Hadi Hasan1, Noel Maalouf2, and Imad H.
Elhajj1 1Authors are with the American University of Beirut, Lebanon
{aaa188, mad51<EMAIL_ADDRESS><EMAIL_ADDRESS>is with the
Lebanese American University, Byblos, Lebanon<EMAIL_ADDRESS>
###### Abstract

Robots have become ubiquitous tools in various industries and households,
highlighting the importance of human-robot interaction (HRI). This has
increased the need for easy and accessible communication between humans and
robots. Recent research has focused on the intersection of virtual assistant
technology, such as Amazon’s Alexa, with robots and its effect on HRI. This
paper presents the Virtual Assistant, Human, and Robots in the loop (VAHR)
system, which utilizes bidirectional communication to control multiple robots
through Alexa. VAHR’s performance was evaluated through a human-subjects
experiment, comparing objective and subjective metrics of traditional keyboard
and mouse interfaces to VAHR. The results showed that VAHR required 41% less
Robot Attention Demand and ensured 91% more Fan-out time compared to the
standard method. Additionally, VAHR led to a 62.5% improvement in multi-
tasking, highlighting the potential for efficient human-robot interaction in
physically- and mentally-demanding scenarios. However, subjective metrics
revealed a need for human operators to build confidence and trust with this
new method of operation.

##  I INTRODUCTION

Virtual assistants have become ubiquitous in our daily lives, integrated into
smartphones, in-car infotainment systems, and smart speakers such as Google
Home and Amazon Alexa. At the same time, robots are increasingly used for
various tasks, from vacuuming homes [1] to delivering goods in hospitals [2].
The potential benefits of integrating virtual assistants and robots are thus
increasingly appealing, as virtual assistants can enhance robots’ physical
capabilities, while robots can provide a natural Human-Robot Interface (HRI)
for virtual assistants. In this paper, we propose VAHR (Virtual Assistant,
Human, and Robots in the loop), a system that connects virtual assistants to
robots through different communication schemes, either directly or indirectly.
Our implementation uses Amazon’s Echo Dot and Alexa virtual assistant, given
their widespread availability and development environment (Alexa Skills Kit
and AWS), but our proposed framework can be easily adapted to other virtual
assistants. We review related work in Section II and describe our system
design in Section III. We then detail the implementation and experimental
setup in Section IV before presenting the results and analysis in Section V.
The conclusion and future work are covered in Section VI.

##  II Related Work

Recent research has focused on developing natural speech interfaces for human-
robot communication using affordable home automation tools like Amazon Alexa.
In this section, we will present various approaches to human-robot voice
interfaces found in the literature, followed by an explanation of the commonly
used HRI standard metrics for evaluating our proposed experiment.

###  II-A Voice-controlled Robots

Earlier attempts of controlling robots through voice date back to 1998 and
2000 in which robots assisted mainly in surgery while responding to a limited
set of voice commands and requiring the surgeon to wear a headset [3, 4].

Recently, robotic assistants such as Lio were deployed in hospitals to assist
in healthcare, and disinfection [5]. Because of Lio’s computational power, its
voice interactions were purely synthesized natively without requiring any
cloud-based solution [5]. However, these native solutions are considerably
more expensive than cloud-based ones. Hence, controlling robots via a friendly
communication scheme has been shifting towards low-cost cloud solutions. This
particular field has attracted many applications in both industry and
academia. In [6], researchers developed EchoBot, an interface between YuMi,
the industrial robot, and Amazon Alexa to facilitate data collection for
Learning from Demonstration (LfD) tasks. To test the efficacy of the
experiment, a user is tasked with guiding YuMi to solve the Tower of Hanoi
puzzle while occupying both hands. The user utters voice commands to the Echo
Dot device to record the state of the two grippers involved in carrying out
the task (open or closed) [6]. This interface was juxtaposed with a keyboard
interface, and the results demonstrated that EchoBot was significantly more
efficient than the keyboard interface [6]. Alexa was also used as a speech
interface in [7] to control the Kuka robotic arm with a gripper. The arm would
respond to voice commands issued by the user to Amazon Alexa, and Alexa relays
the command to the robot to pick and place a wrench from a toolbox via
subscription to an MQTT topic [7].

Further, in [8], a 3D-printed humanoid torso was manipulated through Alexa to
enhance trust in HRI. Additionally, as detailed in [9], a six-degrees-of-
freedom robotic arm was voice controlled through Alexa via a locally-hosted
web server on Ngrok. Ngrok allows a user to host a local server with minimal
effort. In [10], the motors of a semi-autonomous lawn mower were strictly
controlled through the Amazon Tap speaker, a voice automation tool equipped
with Alexa.

A more abstract work, such as in [11], depicts multiple approaches and
architectures to interface one’s robot with Alexa. To our knowledge, no work
in the literature includes an architecture that allows for multi-robot
bidirectional communication and is assessed according to standard HRI metrics.
Our proposed VAHR system is assessed according to the metrics in the following
subsection.

###  II-B Popular HRI Metrics

A comprehensive literature survey of existing HRI metrics over the past two
decades was conducted by Murphy et al. [12]. It can be argued that the most
standardized and cited metrics are summarized in [13]. Saleh et al. [14]
improve these metrics. We list the metrics that are relevant to our study in
what follows.

  * •

Robot Attention Demand (RAD): this metric measures the fraction of time spent
by the user interacting with the robot [14]. Typically, This metric is
calculated by the following formula:

$RAD=\frac{IE}{IE+NT}$ (1)

where NT represents the neglect tolerance, a metric that captures the time
elapsed after the human operator issues the robot command until there is a
notable drop in robot performance. IE denotes the interaction effort
corresponding to the required time to interact with the robot [14]. Both of
the metrics are acquired experimentally.

  * •

Fan-out (FO) evaluates the efficiency and simultaneity of controlling multiple
robots. This usually depicts the operator’s efficiency in operating multiple
homogeneous robots. It is computed by dividing the total task time by the RAD
for a single robot.

$FO=\frac{TotalTaskTime}{RAD}$ (2)

  * •

Trust: to quantify the trust established between a human and a robot, a fuzzy
temporal model is incorporated in [14]. This model takes quantifiable inputs
such as fault size, productivity, and awareness and provides five states of
trust as output (Very Low, Low, Medium, High, and Very High). Each of these
states is associated with a value. (Very Low $\longrightarrow$ 0.1, Low
$\longrightarrow$ 0.3, Medium $\longrightarrow$ 0.5, High $\longrightarrow$
0.7, and Very High $\longrightarrow$ 0.9). The claim in [14] is that trust can
significantly affect RAD. For instance, a human who is less trustworthy of his
or her robot is likely to give greater attention to it even after sending the
command for it to perform independent work. Hence, the following formula comes
about according to [14]

$RAD=DIT+IIT=DIT+NT*(1-T_{r})$ (3)

where DIT is the direct interaction time which is essentially the formula of
(1), indirect interaction time (IIT) is a result of trust (Tr) and neglect
tolerance (NT). In this way, (3) captures the degree of trust the operator
gives to the commanded robot. In controlled experiments, the value of trust
can also be obtained through questionnaires rather than relying on a complex
fuzzy model.

  * •

Degree of mental computation: some tasks may require the operator to
significantly stress his/her short and long-term memory. Steinfeld et al. [13]
mention examples of mental labeling of objects, object-referent association in
working memory, and teleportation tasks. One may resolve this by resorting to
feedback from the robot itself. One of the aims of this project is to allow
robots to report feedback through Amazon Alexa. This metric is also useful in
scenarios to assess controlled experiments. A popular tool named NASA-TLX
(task load index) is typically utilized to measure this metric in controlled
human experiments. It targets six fields in mental and physical workload [15].

  * •

Efficiency and Effectiveness: this is a common metric used in HRI [13].
Efficiency measures the time needed to complete a task, while effectiveness is
the success rate. It calculates the percentage of the robot’s mission that was
successful.

##  III System Design

###  III-A Alexa to Robot Interaction

Having surveyed the different approaches used in the literature, we will
discuss these architectures in light of our approach to the problem.

####  III-A1 MQTT Commands with Shadow Feedback

Figure 1: Our proposed communication scheme utilizes MQTT commands and the
device shadow service.

In this communication scheme, we combine the work done in [16] with the work
done in [17] while abstracting the MQTT broker functionality to AWS. We thus
utilize AWS IoT Core’s MQTT broker to issue commands to robots while providing
feedback to users through Alexa via the device shadow service, as shown in
Fig. 1. Relying only on device shadows as in [17] incurs redundant shadow
requests that increase the cost of the system. Instead, we minimize the number
of shadow requests by using MQTT messages to enable asynchronous messages with
one-to-many communication capability. This ensures efficient scaling of the
communication architecture as more robots are added to the system.

####  III-A2 One-to-Many Communication

With MQTT messages, we can simultaneously broadcast a message from Alexa to
several robots and devices using the publish-subscribe model. By subscribing
to a common MQTT topic, multiple robots can receive commands in parallel as
Alexa publishes these commands. This offloads communication from Lambda to the
MQTT broker and thus ensures timely responses from Lambda. This also allows us
to scale to swarms of robots if needed.

####  III-A3 Asynchronous Communication

Instead of requiring robots to actively query their shadows for commands, we
initiate an MQTT subscriber, which listens to MQTT messages received at a
certain port of the robot. This frees the robot from wasting computations and
bandwidth to poll its device shadow continuously, especially when most
requests are redundant. Now, the robot would only process requests when a
message is received.

We can now describe the communication flow as follows.

  1. 1.

User invokes desired robot-related Alexa skill.

  2. 2.

Desired Alexa Skill interprets the spoken utterance and sends the interpreted
intent to the intent handler endpoint.

  3. 3.

The AWS Lambda function receives the intent as an Alexa Skill intent handler.

  4. 4.

If the invoked intent requires a robot command to be initiated, AWS Lambda
publishes a message with the command values to the respective MQTT topic
(i.e., spin topic to spin a robot).

  5. 5.

All robots subscribed to that topic receive the message with its value when
Lambda publishes to that topic.

  6. 6.

After going through the intent logic and issuing the requested MQTT message
(if needed), the intent handler requests the robot’s shadow to receive the
reported states by the robot (i.e., the command received successfully).

  7. 7.

The handler finally responds with a speech string to be uttered to the user by
the Alexa Skill model, including any updates of the robot’s states, if they
exist.

  8. 8.

On the robot side, the device shadow is updated whenever any of the robot’s
states have changed (i.e., Task Completed or Stuck).

Accordingly, we have created a scalable architecture that allows robots to
integrate robustly and efficiently with Alexa. The advantages and
disadvantages of this service can be summarized as follows.

  * •

Advantages

    1. 1.

Efficient Communication: Fewer shadow requests and MQTT messages are required.

    2. 2.

Flexibility: Enables one-to-many and one-to-one communication between Alexa
and Robots.

    3. 3.

Robustness: The asynchronous nature of MQTT messaging, coupled with the
robustness of the AWS MQTT broker, ensures that the system continues to
function if any of the involved endpoints stops working. Commands and states
would be delayed, but the communication is not fatally broken.

  * •

Disadvantages

    1. 1.

Complexity: Two services must be configured (MQTT Broker and device shadows)
for this setup.

###  III-B Robot to Alexa Voice Interaction

To further validate the importance of using Alexa in robotics, a system is
devised that grants a robot the ability to communicate with Alexa through
speech. This system is constrained to the purposes of the experiment explained
in the next section, but it can be further extended to be deployed in any
application. In essence, the robot utters commands that trigger Alexa to
respond, and according to this response, the robot performs a designated
action sequence.

Figure 2: Robot to Alexa voice System architecture as applied to a weather
application.

To implement this, a collection of tools is utilized. First, Google’s Text to
Speech (gTTS)-API converts assigned text phrases to the robot. For instance,
to call Alexa, the robot would utter a phrase such as ”Alexa, what is the
weather today?” to inquire about the weather. The synthesized speech generated
by gTTS is played via a media player python library dubbed playsound. Alexa
would naturally respond by stating the weather. For the robot to capture
Alexa’s voice utterance, it employs the Speech Recognition engine, which
translates the voice expression to text through Google’s Speech-to-Text API.
This captured text phrase is forwarded to an Amazon Lex chatbot for
interpretation. Amazon Lex is an easy-to-use conversational bot development
tool, part of the AWS suite, capable of recognizing the intent of the text.
According to this intent, the robot would utter back a statement to Alexa
acknowledging its response and executes a certain action accordingly. Figure 2
illustrates the weather example we applied in the experiment if the weather is
sunny. Upon acknowledgment of sunny weather, the robot navigates using the ROS
move_base package to a predefined set of coordinates; a ”sunny” zone.

##  IV Implementation and Experimental Setup

We implemented our proposed design using Python and the AWS IoT SDK. The
system was tested using an Amazon Alexa Echo dot and a mobile robot (placebot)
running ROS melodic. A test setup was created to compare our approach to
traditional robot control methods like a keyboard and mouse. A schematic
diagram describing the testing environment is provided in Fig. 3.

Figure 3: Experiment setup schematic

The experiment involves one test user controlling two robots and finishing
several tasks along the way. The robots (placebots) are mobile robots
programmed to be capable of autonomous navigation inside the allowed
boundaries (white zones). The different robotic arms around the testing
environment are also programmed to be able to pick and place color-coded
(green) items from atop these mobile robots to the ground or from a conveyor
belt (loading zone) to atop one placebot. These functionalities are essential,
as the experiment’s goal is to test our communication architecture and not
robots’ navigation or robot arm placement. Additionally, the testing
environment is primarily virtual, with the robot environments simulated using
ROS and Gazebo. This ensures an ideal and repeatable testing environment,
especially for robot-to-Alexa communication, where a robot has to be near
Alexa to ask and listen to responses. The robots are simulated on a computer
on the same desk where Alexa is placed. The Gazebo world used is shown in Fig.
4.

Figure 4: Gazebo world used for the experiment Figure 5: Experiment setup

The user needs only to order the robots with different commands and solve a
physical jigsaw puzzle. The complete setup is illustrated in Fig. 5, showing
the simulator and the location of the Amazon Echo Dot and the jigsaw puzzle,
all within the user’s reach on the bench. A total of 10 subjects (two females
and eight males) participated in the experiment. Participants’ ages ranged
from 18 to 30 years. The flow of the experiment is explained below.

  1. 1.

The user will go through two runs for each control method, including one
5-minute warm-up run for the user to get used to the controls.

  2. 2.

A test run is deemed complete once all tasks are completed.

  3. 3.

The user is tasked with four tasks: three robot-related tasks and one
distracting task.

  4. 4.

The distracting task entails the user finishing a jigsaw puzzle. No time limit
is set on the tasks.

  5. 5.

The three robot tasks are designed based on functionality:

Task I - Alexa-to-Robot Communication: In this task, a user orders the robots
to move into the loading zone to be loaded with a package (green box), before
ordering them to deliver the package to one of the Zones (A through D). The
zone selection is randomized at the start of the experiment and is briefed to
the user during the experiment (via a display on the computer). This task is
deemed successful if each robot delivers its package to its set Zone. A
failure is marked whenever a robot delivers a package to the wrong room. An
intent is programmed in Alexa to enable the user to order a robot to navigate
to one Zone (Zone A through D or Loading Zone). Robotic arms will load or
unload the robots once they reach the zone waypoint.

Task II - Robot-to-Alexa Communication: In this task, a user orders the robots
to deliver a package based on the weather. The weather is pre-programmed to be
randomized by the Alexa skill and is specific to the experiment (i.e., not the
real weather). As such, the robots are expected to deliver packages to Zone A
for sunny weather and Zone C for rainy weather. There are two ways of
approaching this problem:

     * •

Traditional Control (keyboard and mouse): the user has to ask Alexa for the
weather twice (once for each robot) before ordering navigation to the correct
zones.

     * •

Proposed Communication Architecture: the user can trigger a weather-navigation
intent on Alexa. This will communicate to one robot that the user wants it to
navigate based on the weather. Upon receiving the command, the robot will ask
Alexa for the weather, interpret the response, and autonomously load,
navigate, and unload a package to the correct zone. This will also test the
robot’s ability to communicate with Alexa through voice commands.

This task is deemed completely successful if each robot delivers one package
to the correct zone (based on the weather). Failure in this task is defined as
a robot navigating to the wrong room for the weather (i.e., Zone C when the
weather is sunny). This is marked as a robot’s failure to interpret Alexa’s
response correctly. An intent is programmed in Alexa to enable the user to
order a robot to navigate based on the weather. Robotic arms will load or
unload the robots once they reach the respective zone way-point.

Task III - Robot-to-Robot Communication: In this task, the user is ordered to
make the robots deliver two packages to Zone D in sequential order: placebot 1
delivers a package to Zone D, then placebot 2 delivers another package to zone
D. This task is deemed successful if each robot delivers one package to Zone D
and placebot 1 delivers the package before placebot 2. If the task is
successful, a green LED is turned on inside the simulation environment. The
task is deemed failed if either of the robots fails to deliver the package to
Zone D (i.e., delivered to a different Zone) or the order of delivering the
packages is wrong (i.e., placebot 2 delivers the package before placebot 1).
An intent is programmed in Alexa to enable the user to order a robot to start
communicating with the other robot to coordinate sequential delivery. Robotic
arms will load or unload the robots once they reach the respective zone way-
point.

  6. 6.

Once the experiment is over (all four runs are finished), the user is given a
NASA task load index (TLX) appended with two additional questions specific to
our experiment.

The evaluation criteria used during the experiment are split between
subjective and objective metrics. The subjective metrics are filled in by the
user after they are done with the experiment. The objective metrics are
automatically collected and measured during the experiment. An observer only
fills in the calculated results after each run. The metrics used during the
experiment are summarized in Table I.

TABLE I: Summary of Evaluation Metrics Metric Type | Parameter
---|---
Subjective | Mental Demand
Physical Demand
Temporal Demand
Performance
Effort
Frustration
Confidence
Trust
Objective | Command Success Rate
Task Success Rate
Communication Success Rate
Command Failure Rate
Task Failure Rate
Communication Failure Rate
Robot Attention Demand (RAD)
Fan-out (FO)
Time Taken

##  V Results and Analysis

Objective measures were used to evaluate the ability to communicate with and
command the robots while solving the puzzle. Both methods achieved 100%
command, task, and communication success rates across all trials and for all
test subjects. Four objective metrics (RAD, FO, time taken, and the number of
solved puzzle pieces) were used to quantitatively compare the effect of both
control methods on the users’ performance while completing the assigned tasks.

The results are summarized in Table II per method showing the average score
across both trials. The VAHR method required nearly 41% less Robot Attention
Demand while ensuring almost 93% more Fan-out time. This shows significant
improvement in controlling both robots simultaneously without the need for
increased attention from the user. Although users needed more time to complete
the task using VAHR (254 seconds) compared to the standard mouse-and-keyboard
method (234 seconds), VAHR allowed the users to solve more pieces of the
puzzle averaging a 62.5% improvement. An ANOVA test was conducted to assess
the statistical significance of the obtained objective results. The p-values
obtained from comparing all metrics under both control methods are listed in
Table II. The ANOVA test results show strong evidence against the null
hypothesis (i.e., $p-values<\alpha=0.05$), thus proving that the results are
statistically significant.

TABLE II: Objective Results | VAHR | Mouse & Keyboard | p-value
---|---|---|---
Interaction Time (s) | 111 | 175 | 2.9e-4
Neglect Time (s) | 142 | 59 | 2e-5
RAD [0, 1] | 0.44 | 0.75 | 7.3e-7
Fan-out (s) | 578 | 314 | 6.3e-4
Solved Puzzle Parts | 13 | 8 | 3e-2

Figure 6 shows the averages of the NASA Task Load Index subscales. The NASA
Task Load Index consists of a scale of 21 assessment levels which we converted
to a percentile score from 0% at level 1 to 100% at level 21. The lower score
reflects less load on the subject under testing. While VAHR required 17.6%
more mental demand on average than the traditional communication method, the
physical load on the subject was decreased by 41.6%. The temporal load was
almost the same for both methods with a 1% difference. The users assess their
own performance to be 22.6% better when using the traditional method which is
also reflected in the 9.6% more trust. The subjective assessment of
performance contradicts the objective evaluation scores as the subjects show
superior performance while completing the tasks using VAHR. We estimate that
the discrepancy goes back to the fact that the proposed VAHR method is
relatively new to the subjects who need time to build trust in the new control
method. VAHR required 3.5% more effort on average from the users to complete
the assigned tasks but decreased the frustration level by 5.1% and improved
the subjects’ confidence in the ability to complete the tasks successfully by
1.7%. The ANOVA test conducted on the subjective results obtained from both
control methods resulted in a p-value of 0.28 ($p-value<\alpha=0.05$) which
doesn’t show sufficient statistical significance for the obtained results.
Individual ANOVA tests were conducted on the results of each subjective metric
of the NASA TLX to check for statistical significance. None of the subjective
results received a p-value lower than 0.05 across the two methods and
therefore the null hypothesis cannot be rejected for any of the subjective
metrics.

Figure 6: Subjective results comparison

Since VAHR is relatively new to the users, we consider the change in the NASA
Task Load Index scores across different trials. As shown in Fig. 7, the users
felt that the task required less mental and temporal demand in the second
trial compared to the first one, with each of the two metrics registering an
improvement of around 29%. Improvement is also evident in the effort index,
with around 15% and the highest percentage improvement is in the frustration
index which reaches more than 40% in the second trial. On the other hand, VAHR
required 29% more physical demand on average from the users in their second
trial where they also felt that their performance and confidence dropped by
around 3% and 8% respectively. The trust score remained constant across both
trials.

Figure 7: VAHR - subjective results improvement across trials

##  VI Conclusion

We present VAHR, a scalable cloud-based communication paradigm that connects
humans to robots through virtual assistants. This communication system enables
all agents to communicate freely with each other through speech. While it
requires more setup than traditional approaches, cloud services provide
increased robustness and scalability. We test our system against standard
mouse-and-keyboard control through a set of human-subject experiments, which
demonstrated the advantage of using VAHR in objective metrics, providing at
least 41% improvement. Part of the subjective measures noted by users shows a
preference for mouse-and-keyboard control, which we attribute to unfamiliarity
with this novel communication scheme. However, the high p-values obtained from
the ANOVA tests on both combined and individual subjective metrics show low
statistical significance of the subjective test results. This adds evidence to
the unfamiliarity of the subjects with the proposed control/communication
method, even from a self-assessment point of view. We suggest that our system
would be best suited for environments such as warehouses and hospitals, where
physical tasks have to be balanced with robot control.

## References

  * [1] “Roomba® robot vacuums.” [Online]. Available: https://www.irobot.com/en˙US/roomba.html
  * [2] “Tug autonomous mobile robots for healthcare and hospitality.” [Online]. Available: https://aethon.com/products/
  * [3] S. Okada, Y. Tanaba, H. Sugawara, H. Yamauchi, S. Ishimori, and S. Satoh, “Thoracoscopic major lung resection for primary lung cancer by a single surgeon with a voice-controlled robot and an instrument retraction system,” _The Journal of Thoracic and Cardiovascular Surgery_ , vol. 120, no. 2, p. 414–415, 2000.
  * [4] L. Mettler, M. Ibrahim, and W. Jonat, “One year of experience working with the aid of a robotic assistant (the voice-controlled optic holder aesop) in gynaecological endoscopic surgery,” _Human Reproduction_ , vol. 13, no. 10, p. 2748–2750, 1998.
  * [5] J. Mišeikis, P. Caroni, P. Duchamp, A. Gasser, R. Marko, N. Mišeikienė, F. Zwilling, C. de Castelbajac, L. Eicher, M. Früh, and H. Früh, “Lio-a personal robot assistant for human-robot interaction and care applications,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 5339–5346, 2020.
  * [6] R. Kapadia, S. Staszak, L. Jian, and K. Goldberg, “Echobot: Facilitating data collection for robot learning with the amazon echo,” _2017 13th IEEE Conference on Automation Science and Engineering (CASE)_ , 2017.
  * [7] J. Jungbluth, K. Siedentopp, R. Krieger, W. Gerke, and P. Plapper, “Combining virtual and robot assistants a case study about integrating amazon’s alexa as a voice interface in robotics,” _University of Luxembourg_.
  * [8] B. A. Erol, C. Wallace, P. Benavidez, and M. Jamshidi, “Voice activation and control to improve human robot interactions with iot perspectives,” in _2018 World Automation Congress (WAC)_ , 2018, pp. 1–5.
  * [9] Z. Wang, D. Chen, and P. Xiao, “Design of a voice control 6dof grasping robotic arm based on ultrasonic sensor, computer vision and alexa voice assistance,” in _2019 10th International Conference on Information Technology in Medicine and Education (ITME)_ , 2019, pp. 649–654.
  * [10] J. A. Solorio, J. M. Garcia-Bravo, and B. A. Newell, “Voice activated semi-autonomous vehicle using off the shelf home automation hardware,” _IEEE Internet of Things Journal_ , vol. 5, no. 6, pp. 5046–5054, 2018.
  * [11] A. Hidalgo-Paniagua, A. Millan-Alcaide, J. P. Bandera, and A. Bandera, “Integration of the alexa assistant as a voice interface for robotics platforms,” in _Robot 2019: Fourth Iberian Robotics Conference_ , M. F. Silva, J. Luís Lima, L. P. Reis, A. Sanfeliu, and D. Tardioli, Eds. Cham: Springer International Publishing, 2020, pp. 575–586.
  * [12] R. Murphy and D. Schreckenghost, “Survey of metrics for human-robot interaction,” in _Proceedings of the 8th ACM/IEEE International Conference on Human-Robot Interaction_ , ser. HRI ’13. IEEE Press, 2013, p. 197–198.
  * [13] A. Steinfeld, T. Fong, D. Kaber, M. Lewis, J. Scholtz, A. Schultz, and M. Goodrich, “Common metrics for human-robot interaction,” in _Proceedings of the 1st ACM SIGCHI/SIGART Conference on Human-Robot Interaction_ , ser. HRI ’06. New York, NY, USA: Association for Computing Machinery, 2006, p. 33–40. [Online]. Available: https://doi.org/10.1145/1121241.1121249
  * [14] J. A. Saleh and F. Karray, “Towards generalized performance metrics for human-robot interaction,” in _2010 International Conference on Autonomous and Intelligent Systems, AIS 2010_ , 2010, pp. 1–6.
  * [15] S. G. Hart and L. E. Staveland, “Development of nasa-tlx (task load index): Results of empirical and theoretical research,” _Human mental workload_ , vol. 1, no. 3, pp. 139–183, 1988.
  * [16] A. Hidalgo-Paniagua, A. Millan-Alcaide, J. P. Bandera, and A. Bandera, “Integration of the alexa assistant as a voice interface for robotics platforms,” in _Iberian Robotics conference_. Springer, 2019, pp. 575–586.
  * [17] J. Jungbluth, R. Krieger, W. Gerke, and P. Plapper, “Combining virtual and robot assistants-a case study about integrating amazon’s alexa as a voice interface in robotics,” in _Robotix-Academy Conference for Industrial Robotics (RACIR) 2018_. Shaker, 2018, p. 5.


# Modelling thermoresponsive polymer brush by mesoscale computer simulations

D. Yaremchuk  $|$ , O. Kalyuzhnyi  $|$ , J. Ilnytskyi  $|$ Corresponding
author<EMAIL_ADDRESS>
(Received April 17, 2023)

###### Abstract

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.

Key words: thermoresponsive polymers, PNIPAM, dissipative particle dynamics

###### Abstract

Ми розглядаємо функціональну поверхню, що містить термочутливі полімерні
ланцюги, матеріал, який широко застосовується в технології та біомедицині. З
метою опису її характерних просторово-часових масштабів, комп’ютерне
моделювання виконується із залученням огрублених мезоскопічних моделей.
Використана тут модель була запропонована авторами роботи [Soto-Figueroa та
ін., Soft Matter, 8, 1871 (2012)], і вона описує основну особливість
полі(N-ізопропілакриламіду) (PNIPAM), а саме – швидку зміну його
гідрофільності при переході через нижню критичну температуру розчинності
(LCST). Для випадку ізольованого ланцюга ми досліджуємо властивості
масштабування його радіусу гірації, відстані між першим та останнім
мономерами, низку функцій розподілу вказаних характеристик та профілі густини
мономерів залежно від відстані від поверхні як нижче, так і вище за LCST. Для
випадку моделі термочутливої полімерної щітки дослідження концентруються на
пошуку оптимальної густини пришпилення, за якої зміна висоти щітки при
переході через LCST досягає максимального значення. Надано інтерпретацію
механізму зміни конформацій ланцюжків щітки в термінах блобів Александера-де
Жена та ступеня сольватації полімерних ланцюгів у щітці.

Ключов слова: термочутливі полімери, PNIPAM, дисипативна динаміка

##  1 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^{\circ}$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. This indicates high sensitivity of the coil-to-globule transition to
the fine details of the force fields. Besides the obvious role of the
temperature, the pressure also plays an important role in this transition.
This effect was addressed by Tavagnacco et al. [14]. In particular, a coil-to-
globule transition was found to occur up to large pressures, and a new high
pressure globular state was found. It is characterized by a more structured
hydration shell that is closer to PNIPAM hydrophobic domains, as compared to
the globular state at atmospheric pressure.

A single PNIPAM chain, grafted inside a cylindrical pore that resembles a
carbon nanotube, was considered by Alaghemandi et al. [15]. To suppress the
collapse of PNIPAM on a substrate, they reduced the AMBER force field
interaction strengths by a factor of ten. For this case, the end-to-end
distances, radii of gyration, density profiles, the number of hydrogen bonds,
and radial distribution functions for a PNIPAM chain were evaluated, in order
to characterize the coil-to-globule transition on a level of the monomers
rearrangement. The cases of the PNIPAM containing a surfactant molecule [16]
and PNIPAM based hydrogels [17, 18] have also been studied. In particular,
Walter et al. [17] reported a fair agreement between the experimental and
computer simulation studies of the coil-to-globule transition in the PNIPAM
hydrogel in water.

Atomic resolution computer simulations provide invaluable information on the
chemistry detailed properties of the PNIPAM chains and underlying mechanisms
of their conformational changes, as well as their hydration and specific
features of dehydration. However, most applications of the PNIPAM based
systems occur at a relatively large length and time scales, be these either
micellar systems for temperature controlled drug delivery or polymer brushes
for smart commanding surfaces. Atomistic simulations at this scale in most
cases prove to be too costly, which calls for more economical coarse-grained
modelling approaches. These treat each characteristic group of atoms as a
single particle, and their interaction potentials represent total interaction
of their constituent atoms with their surroundings. This enables a drastic
reduction of the degrees of freedom, but should be appropriately parametrized
to retain the principal features of the initial PNIPAM chain. General concepts
of a rigorous coarse-graining of polymers is covered in a number of reviews,
e.g., in the recent one by Dhamankar and Webb [19].

Mid-level coarse-graining of the PNIPAM polymer involves treating its
repeating unit as a trimer. It includes the hydrocarbon group of the main
chain and a hydrophilic amide group and a hydrophobic isopropyl group of a
side chain (see, e.g., references [20, 21]). Hence, each such group contains
2–3 heavy atoms. Abbott et al. [20] followed the coarse-graining recipe by
Shinoda et al. [22], when the nonbonded interactions were fit to experimental
thermodynamic data. It was necessary to derive different potentials at both
temperatures, $280\,\mathrm{K}$ and $330\,\mathrm{K}$, to capture a correct
temperature-dependence of the interactions. The coil-to-globule transition was
observed for the longer, $N=18$ and $N=30$, oligomers only. Botan et al. [23]
developed an implicit solvent coarse-grained model coupled to a lattice
Boltzmann hydrodynamics, but accounting for these interactions was
insufficient to capture the coil-to-globule transition. Pérez-Ramírez and
Odriozola [21] developed a similar level coarse-grained model based on the
Martini force field. The amide moieties include an electric dipole, and
polymer chains of different sizes accurately reproduce the thermal behaviour
without including temperature-dependent parameters. Bejagam et al. [24]
applied a data-driven machine-learning method combined with a coarse-grained
modelling to study the conformations of PNIPAM in solutions. The temperature
independent model being developed accurately predicts the LCST, at the same
time retaining the tacticity in the presence of an explicit water model. The
nonmetric multidimensional scaling method shows the presence of multiple
metastable states of PNIPAM during its coil-to-globule transition above the
LCST. All in all, the coarse-grained modelling of this level has its strong
points in drastic reduction of the degrees of freedom retaining the molecular
topology of PNIPAM polymer and relevant details of its chemical moieties, and,
therefore, may be seen as an optimum compromise in modelling PNIPAM as for
now.

It is, however, tempting, to go one step further and to examine the
capabilities of yet more coarse-grained approach, like the one used in
references [25, 26]. In this case, PNIPAM is represented as a linear chain.
The logic for employing this approach is at least twofold. Firstly, the
experimentally observed effect may occur on a rather large length scale and
its purely physical, rather than chemical, interpretation might be sufficient.
This is the case, for instance, for the setup considered in a recent paper by
Minko group [27]. Here, the separation of colloid particles from a solution is
achieved by a local oscillating repulsive force generated at the
microstructured stimuli-responsive polymer interface, in order to switch
between adsorption and mechanical-force-facilitated desorption of the
particles. The polymer interface is 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 $\chi$
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 $\chi$ 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.

##  2 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 $k_{\rm B}T=1$, where $k_{\rm B}$ is the Boltzmann constant
and $T$ is the temperature, time unit is $t^{*}=1$. The monomers are bonded
via harmonic springs which results in the bonding force:

$\mathbf{F}^{B}_{ij}=-k\mathbf{x}_{ij},$ (1)

where $\mathbf{x}_{ij}=\mathbf{x}_{i}-\mathbf{x}_{j}$, $\mathbf{x}_{i}$ and
$\mathbf{x}_{j}$ are the positions of the $i$-th and $j$-th bead, and $k$ is
the associated energy. The total non-bonded force $\mathbf{F}_{ij}$ acting on
$i$-th bead from its $j$-th counterpart can be expressed as a sum of three
contributions:

$\mathbf{F}_{ij}=\mathbf{F}^{\mathrm{C}}_{ij}+\mathbf{F}^{\mathrm{D}}_{ij}+\mathbf{F}^{\mathrm{R}}_{ij},$
(2)

where $\mathbf{F}^{\mathrm{C}}_{ij}$ is the conservative force caused by the
repulsion of the $i$-th and $j$-th beads, $\mathbf{F}^{\mathrm{D}}_{ij}$ is
the dissipative force responsible for friction between beads (similar to the
Stokes force), and the random force $\mathbf{F}^{\mathrm{R}}_{ij}$ which,
together with the dissipative force, maintains a constant temperature of the
system. Expressions for these terms are pretty standard and are given below:

$\mathbf{F}^{\mathrm{C}}_{ij}=\left\\{\begin{array}[]{ll}a(1-x_{ij})\displaystyle\frac{\mathbf{x}_{ij}}{x_{ij}},&x_{ij}<1,\\\
0,&x_{ij}\geqslant 1,\end{array}\right.$ (3)
$\mathbf{F}^{\mathrm{D}}_{ij}=-\gamma
w^{\mathrm{D}}(x_{ij})(\mathbf{x}_{ij}\cdot\mathbf{v}_{ij})\frac{\mathbf{x}_{ij}}{x^{2}_{ij}},$
(4) $\mathbf{F}^{\mathrm{R}}_{ij}=\sigma
w^{\mathrm{R}}(x_{ij})\theta_{ij}\Delta
t^{-1/2}\frac{\mathbf{x}_{ij}}{x_{ij}},$ (5)

where $x_{ij}=|\mathbf{x}_{ij}|$,
$\mathbf{v_{ij}=\mathbf{v_{i}}}-\mathbf{v_{j}}$, $\mathbf{v_{i}}$ is the
velocity of the $i$-th monomer, and $a$ is the amplitude for the conservative
repulsive force. The dissipative force has an amplitude $\gamma$ and its
decrease with the particle separation is given by the function
$w^{\mathrm{D}}(x_{ij})$. Similarly, the amplitude for the random force is
$\sigma$ and its decrease is provided by $w^{\mathrm{R}}(x_{ij})$.
$\theta_{ij}$ is a Gaussian random variable. As shown in [33], to satisfy the
requirement of a detailed balance, the relation $\sigma^{2}=2\gamma$ and
$w^{\mathrm{D}}(x_{ij})=\large[w^{\mathrm{R}}(x_{ij})\large]^{2}$ should hold.
Here, we use a quadratic decaying form for the weight function:

$w^{\mathrm{D}}(x_{ij})=\left[w^{\mathrm{R}}(x_{ij})\right]^{2}=\left\\{\begin{array}[]{ll}(1-x_{ij})^{2},&x_{ij}<1,\\\
0,&x_{ij}\geqslant 1.\end{array}\right.$ (6)

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 $N$ such monomers that are bonded pairwise via the force (1) with the
spring constant $k=4$.

Figure 1: (Colour online) Schematic representation of the PNIPAM chain
fragment modelled as a linear chain of coarse-grained beads.

In general, repulsive strength of the conservative force, $a$, 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 $a$ was found from matching the inverse
compressibility of a model to that of water at normal conditions [32]. This
resulted in the value of $a=25$. We denote polymer beads by a subscript “$p$”,
and those representing water, by “$w$”. It is assumed that for the polymer-
polymer and water-water interaction, one has $a_{pp}=a_{ww}=25$. Mixed,
polymer-water interaction strength, $a_{pw}$, depends on the polymer level of
hydrophilicity. In the case of PNIPAM, one expects $a_{pw}\approx 25$ below
LCST and $a_{pw}>25$ above it. Hence, the strength of the polymer-water
hydrogen bonding does enter this model implicitly, via the value of the
$a_{pw}$. The temperature dependence of $a_{pw}$ was estimated from that for
the Flory-Huggins mixing parameter $\chi$, 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 $T=298\,\textrm{K}$ (below LCST), the estimate for $a_{pw}$
is $a_{pw}\approx 25.6$, whereas at $T=310\,\textrm{K}$ (above LCST), one has
$a_{pw}\approx 38$.

The simulation setup is as follows. The simulation box dimensions are $L_{x}$,
$L_{y}$ and $L_{z}$, 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 $z=0$ and $z=L_{z}$.
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, $\mathbf{x}_{ij}$ denotes the vector that runs from the
center of $i$-th bead towards the nearest wall and is normal to it. The
amplitude of the bead-wall repulsion is fixed at $a_{pw}=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 $M$ of chains, and one of the end
beads of each chain is grafted to a specified fixed point of the wall located
at the $z=0$ plane. Grafting is performed by using the harmonic force (1) with
$k=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 $M$. Specifically, the grid of $K\times K$
grafting points is made in the form of a square lattice, where $K$ is the
integer value of $\sqrt{M}$. The rest $M-K^{2}$ grafting points are
distributed randomly across the wall. The grafting density is measured as the
surface number density of grafts, $\rho_{g}=M/(L_{x}L_{y})$. 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 $\Delta t=0.04$ in the DPD
units. For the case of an isolated chain, $M=1$, considered in section 3,
$8\cdot 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, $M>1$, covered in section 4, $10^{6}$ DPD steps were used, where first
$3\cdot 10^{5}$ of them were discarded for relaxation purpose.

##  3 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 polymer-water interactions specified by the magnitude of $a_{pw}$ in
the expression of a soft repulsive force (3). Below LCST, the value of
$a_{pw}\approx 25.6$ is found [28], which is pretty close to the case of
athermal solvent, $a_{pw}=a_{pp}=a_{ww}=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 $T=298\,\mathrm{K}$ (below
LCST) and at $310\,\mathrm{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 $z=0$ and
surrounded by a solvent that represents water. The details of modelling are
provided in section 2. The set of chain lengths, $N=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:
$L_{x}=L_{y}=L_{z}=b_{0}N^{0.8}$, where $b_{0}\approx 0.9$ [38] is an
equilibrium bond length between the monomers for the chosen values of $k$ and
$a$ 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\cdot 10^{6}$ DPD steps for each chain length $N$, 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,
$\mathbf{G}$, with its components defined as [43]

$G_{\alpha\beta}=\frac{1}{N}\sum_{i=1}^{N}(r_{i,\alpha}-R_{\alpha})(r_{i,\beta}-R_{\beta}).$
(7)

Here, $r_{i,\alpha}$ and $R_{\alpha}$ denote Cartesian coordinates of $i$-th
monomer and these for the centre of mass of a chain, respectively, and
$\alpha=1,2,3$.

Figure 2: (Colour online) Scaling of the average radius of gyration $R_{g}$ of
a single model PNIPAM chain with the number of monomers, $N$, at two
temperatures, $T=298\,\mathrm{K}$ and $310\,\mathrm{K}$. DPD data are
displayed via symbols. Both cases of a free and grafted chains are shown and
the results for their fits to the power law (8) are provided in the plots.

The eigenvalues $\lambda_{\alpha}$ of the gyration tensor are sorted as
follows $\lambda_{1}>\lambda_{2}>\lambda_{3}$. These allow one to evaluate the
instantaneous squared radius of gyration of a chain,
$R_{g}^{2}=\sum_{\alpha=1}^{3}\lambda_{\alpha}$. The average value for the
radius of gyration is defined as follows:

$\bar{R}_{g}=\langle R_{g}^{2}\rangle^{1/2}\sim l\,N^{\nu},$ (8)

where the averaging is performed over the time trajectory and we already used
the scaling law with respect to the number of monomers, $N$, [37, 36]. Here,
$l$ is a pre-factor, and $\nu$ is the scaling exponent which depends on the
solvent quality. In particular, $\nu\approx 0.5882$ [44] for the case of good
solvent, $\nu=1/2$ for the random walk (Gaussian chain) and $\nu=1/3$ for the
case of poor solvent. Grafted chain is expected to obey the power law with the
same exponent $\nu$ as its free counterpart but with a different pre-factor
$l$ [45].

Simulation data obtained for $\bar{R}_{g}$ 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 $\bar{R}_{g}$ is essentially
higher at $T=298\,\mathrm{K}$ than at $310\,\mathrm{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 $\nu$. The estimates for $\nu\approx
0.582-0.586$ at $T=298\,\mathrm{K}$ are in a very good agreement with the
known most accurate result $\nu\approx 0.5882$ [44] for the Flory exponent. At
$T=310\,\mathrm{K}$, however, the value of $\nu\approx 0.300$ is close to but
lower than the expected exact result, $\nu=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 $l$ but the difference is
rather minor, i.e., within the accuracy of simulations.

Figure 3: (Colour online) Scaling of the chain length $R_{e}$ of a single
model PNIPAM chain with the number of monomers, $N$, at two temperatures,
$T=298\,\mathrm{K}$ and $310\,\mathrm{K}$. DPD data are displayed via symbols.
Both cases of a free and grafted chains are shown and the results for their
fits to the power law (9) are provided in the plots.

Similar results are obtained for the average end-to-end distance of a chain
defined via

$\bar{R}_{e}=\langle R_{e}^{2}\rangle^{1/2}\sim l\,N^{\nu},$ (9)

where $R_{e}^{2}=\sum_{\alpha=1}^{3}(r_{1,\alpha}-r_{N,\alpha})^{2}$ and the
averaging is performed over the time trajectory. Here, the pre-factor $l$ and
the critical exponent $\nu$ have the same meaning as in equation (8). In
general, the end-to-end distance is found to be more demanding in terms of the
simulation length and it displays a scaling behaviour starting from longer
chain lengths, $N\geqslant 20$. The estimates for $\nu$ are acceptable for the
case of $T=298\,\mathrm{K}$, whereas the value $\nu\approx 0.27$ at
$T=310\,\mathrm{K}$ is even lower than that obtained from the analysis of
$\bar{R}_{g}$.

Figure 4: (Colour online) Left-hand frame: scaled distribution
$N^{\nu}\phi(R^{*}_{g})$ for the dimensionless radius of gyration
$R^{*}_{g}=R_{g}/\bar{R}_{g}$ at $T=298\,\mathrm{K}$. Here, the respective
scaling laws for $\bar{R}_{g}$, as found in the left-hand frame of figure 2,
are used. Right-hand frame: scaled distributions
$N^{-\alpha}\phi(R^{\prime}_{g})$ for dimensionless radius of gyration
$R^{\prime}_{g}$ at $T=310\,\mathrm{K}$. The latter is obtained by shifting
the $R_{g}$: $R^{\prime}_{g}=1+(R_{g}-\bar{R}_{g})N^{\beta}$. The pair of
unknown exponents $(\alpha,\beta)$ are found numerically as $(0.19,0.3)$ and
$(0.25,0.36)$ for a grafted and free chain, respectively, and the scaling laws
for $\bar{R}_{g}$ are used as shown in the right-hand frame of figure 2.

We switch our attention now to the distributions of the $R_{g}$ values during
the simulation time. For the good solvent case, $T=298\,\mathrm{K}$, such
distribution can be written in terms of the dimensionless radius of gyration,
$R^{*}_{g}=R_{g}/\bar{R}_{g}$, where the scaling law for $\bar{R}_{g}$ is
provided in figure 2. The data for different chain length $N$ obey the same
master curve for the scaled distributions

$N^{\nu}\phi(R^{*}_{g})=f_{L}(R^{*}_{g})=A_{L}\exp\left(-\frac{1}{[R^{*}_{g}]^{\frac{3}{3\nu-1}}}-[R^{*}_{g}]^{\frac{1}{1-\nu}}\right),$
(10)

as shown in figure 4. The function $f_{L}(x)$ is suggested by Lhuillier [46],
following his considerations on the probabilities of fully collapsed and fully
stretched conformations. Here, $A_{L}\approx 42$ is the pre-factor, the first
term within the parenthesis is dominant at small $R^{*}_{g}$, whereas the
second term defines the rate at which $f_{L}(R^{*}_{g})$ decays at large
$R^{*}_{g}$. 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,
$T=310\,\mathrm{K}$. We observed that the initial distributions, in terms of
$R_{g}$, shift towards larger $R_{g}$ and increase moderately in their height
upon the increase of the chain length, $N$. Therefore, the following
expression for the shifted radius of gyration, $R_{g}^{\prime}$, is proposed
in this case (note that in the DPD simulations the length unit is assumed to
be equal to unity)

$R^{\prime}_{g}=1+(R_{g}-\bar{R}_{g})N^{\beta},$ (11)

and the distribution scales as $N^{-\alpha}\phi(R^{*}_{g})$. The pair of
unknown exponents, $(\alpha,\beta)$, is obtained by numeric fitting, yielding
the values $(0.19,0.3)$ and $(0.25,0.36)$ for a grafted and free chain,
respectively. The scaling laws shown in the right-hand frame of figure 2, are
used for $\bar{R}_{g}$. The results for scaled distributions
$N^{-\alpha}\phi(R^{*}_{g})$ are presented in the right-hand frame of figure 4
and one may note a rather good match of the simulation data for different $N$.
The respective master curves at both temperatures, $T=298\,\mathrm{K}$ and
$310\,\mathrm{K}$, are remarkably similar in their shape, but we did not
explore this similarity in any deeper terms here.

Figure 5: (Colour online) Scaled distribution $N^{\nu}\psi(R^{*}_{e})$ at
$T=298\,\mathrm{K}$. Here, $R^{*}_{e}=R_{e}/\bar{R}_{e}$, where the scaling
laws for $\bar{R}_{e}$ are given in the right-hand frame of figure 3.

Now, we discuss similar distributions for the end-to-end distance, $R_{e}$.
Here, only the case of a good solvent, $T=298\,\mathrm{K}$, is considered. The
dimensionless end-to-end distance, $R^{*}_{e}=R_{e}/\bar{R}_{e}$, is
introduced, where the scaling laws for $\bar{R}_{e}$ are shown in the left-
hand frame of figure 3. Similarly to the case of the $R_{g}$, scaling of the
distributions $N^{\nu}\psi(R^{*}_{e})$, obtained at different chain length
$N$, enables us to introduce the master curve

$N^{\nu}\psi(R^{*}_{e})=f_{C}(R^{*}_{e})=A_{C}[R^{*}_{e}]^{\delta}\exp\left(-[s_{C}R^{*}_{e}]^{\frac{1}{1-\nu}}\right),$
(12)

see figure 5. The analytic form for the function $f_{C}(x)$ is constructed
following the analysis by des Cloizeaux [47] and de Gennes [48]. Using numeric
fits for the simulation data, we found that $A_{C}\approx 5.65$,
$\delta\approx 2.63$ and $s_{C}\approx 1.15$. Since $\delta>1$, the shape of
$N^{\nu}\psi(R^{*}_{e})$ at small $R^{*}_{e}$ is convex. One should mention
that the theoretical analysis [47, 48] yields the value of
$\delta=(1-\gamma)/\nu\approx 0.27<1$ resulting in a concave curve at small
$R^{*}_{e}$. We believe that this discrepancy is the result of the excluded
volume effects at small $R_{e}$ in the DPD modelling approach, and these
reduce the probability for the values $R^{*}_{e}\to 0$.

Figure 6: (Colour online) Left-hand frame: scaled density profiles
$N^{\nu}\rho_{z}(z^{*})$ for grafted and free chain, respectively at
$T=298\,\mathrm{K}$. Here, $z^{*}=z/\bar{R}_{g}$, where the scaling law for
$\bar{R}_{g}$ follows that of figure 2. Right-hand frame: the same at
$T=310\,\mathrm{K}$. Symbols represent DPD simulation data, solid lines,
marked $f_{g}(x)$ and $f_{f}(x)$, are the fits to the analytic forms (13) and
(14).

Finally, we examine the density profile $\rho_{z}(z)$ for a single model
PNIPAM chain grafted to the plane $z=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 $T=298\,\mathrm{K}$, the density profile can be built in terms
of the dimensionless separation from the wall, $z^{*}=z/\bar{R}_{g}$ using the
same scaling law for $\bar{R}_{g}$ as above. The data for a scaled profile,
$N^{\nu}\rho_{z}(z^{*})$, obtained at different $N$, 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 $z=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
$N^{\nu}\rho_{z}(z^{*})$ can be constructed following considerations by de
Gennes [37]. In particular, at small $z$, one obtains
$N^{\nu}\rho_{z}(z^{*})\sim(z^{*})^{\frac{1-\nu}{\nu}}$. At large $z^{*}$, we
assume an exponential decay of $\rho_{z}(z^{*})$ similarly to equation (10),
but with unknown parameters. Then, we assume that both terms can be multiplied

$N^{\nu}\rho_{z}(z^{*})=f_{g}(z^{*})=A_{g}(z^{*})^{\frac{1-\nu}{\nu}}\exp{\left[-f_{g}\,(z^{*})^{\alpha_{g}}\right]}.$
(13)

Numeric fits of simulation data yield $A_{g}\approx 0.3$, $f_{g}\approx 0.56$,
and $\alpha_{g}\approx 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 $\rho_{z}(z^{*})$ in this case can be constructed by
modifying the equation (13) by: (i) removing the term
${z^{*}}^{\frac{1-\nu}{\nu}}$, and, (ii) shifting the argument in the
exponential term

$N^{\nu}\rho_{z}(z^{*})=f_{f}(z^{*})=A_{g}\exp{\left[-f_{g}(z^{*}+s_{f})^{\alpha_{f}}\right]}.$
(14)

Note, that the same coefficients, $A_{g}$, $f_{g}$, and the exponent
$\alpha_{g}$ enter this expression, whereas the amount of shift is given by
$s_{f}\approx 0.247$. This analytic form very well describes the scaled
profile $\rho_{z}(z^{*})$ for a free chain, too, see the left-hand frame of
figure 6.

The situation is quite similar at $T=310\,\mathrm{K}$, in which case the
scaling law for $\bar{R}_{g}$, shown in the right-hand frame of figure 2, is
used. The DPD simulation data for both grafted and free chain, represented in
terms of $z^{*}=z/\bar{R}_{g}$, reasonably well follow their respective
universal shapes $\rho_{z}(z^{*})$, see right-hand frame of figure 6. One
should remark that both shapes are very similar to their counterparts at
$T=298\,\mathrm{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 $\bar{R}_{g}$. Similarly to the case of figure 4, we have not performed
analytic fits for the $\rho_{z}(z^{*})$ distributions at $T=310\,\mathrm{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 $a_{pw}$ between PNIPAM monomers and the
beads representing water, see equation 3. Both cases of $T=298\,\mathrm{K}$
(below the LCST), and $T=310\,\mathrm{K}$ (above it) are considered. At the
first stage of this analysis, the scaling laws for the average radius of
gyration, $\bar{R}_{g}$, and these for the end-to-end distance, $\bar{R}_{e}$,
are obtained by performing simulations at a range of chain lengths, $N=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 $R_{g}$ and $R_{e}$ 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, $\bar{R}_{g}$ and $\bar{R}_{e}$, their
distributions, and the density profile of monomers in relation to the wall.

##  4 Thermoresponsive polymer brush

We proceed now to the case of a polymer brush formed of $M$ grafted chains
according to the procedure explained in detail in section 2. The length of all
chains is fixed at $N=60$. The new parameter that is brought into a play here
is the grafting number density of chains, $\rho_{g}$, introduced in section 2.
At low $\rho_{g}$, 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 $\rho_{g}$, 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 $\rho_{g}$. Scaling properties of polymer chains
are widely covered in a number of studies [37, 49, 50, 51, 52, 34, 53, 45].
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]

$\bar{h}=2\frac{\int_{0}^{L_{z}}z\rho_{z}(z)\mathrm{d}z}{\int_{0}^{L_{z}}\rho_{z}(z)\mathrm{d}z},$
(15)

where $\rho_{z}(z)$ is the density profile for the brush beads, and $L_{z}$ 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 $\bar{R}_{g}$ and $\bar{R}_{e}$, are already introduced in equations (8)
and (9). We also complement our analysis by the shape characteristics such as:
the average asphericity, $A$, and the shape descriptor, $B$, [39, 41, 42]
defined as

$\bar{A}=\frac{1}{6}\sum_{i=1}^{N}\left\langle\frac{\sum_{\alpha=1}^{3}(\lambda_{\alpha}-\bar{\lambda})^{2}}{\bar{\lambda}^{2}}\right\rangle,\qquad\bar{B}=3\sum_{i=1}^{N}\left\langle\frac{\bar{\lambda}-\lambda_{2}}{\lambda_{1}}\right\rangle,$
(16)

where $\bar{\lambda}=\frac{1}{3}\sum_{\alpha=1}^{3}\lambda_{\alpha}$, 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 $B=-1$ for an infinitely thin disc, $B=1$ for infinitely
thin rod and $B=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 $T=298\,\textrm{K}$ and $310\,\textrm{K}$, we
consider the set of ratios

$R_{g}^{\dagger}=\frac{\bar{R}_{g}(T=298\,\textrm{K})}{\bar{R}_{g}(T=310\,\textrm{K})},~{}~{}~{}R_{e}^{\dagger}=\frac{\bar{R}_{e}(T=298\,\textrm{K})}{\bar{R}_{e}(T=310\,\textrm{K})},~{}~{}~{}A^{\dagger}=\frac{\bar{A}(T=298\,\textrm{K})}{\bar{A}(T=310\,\textrm{K})},~{}~{}~{}B^{\dagger}=\frac{\bar{B}(T=298\,\textrm{K})}{\bar{B}(T=310\,\textrm{K})},$
(17)

between respective values obtained at these two temperatures.

Figure 7: (Colour online) Characteristic ratios, $h^{\dagger}$,
$R_{g}^{\dagger}$, $R_{e}^{\dagger}$, $A^{\dagger}$, and $B^{\dagger}$ as the
functions of the grafting density $\rho_{g}$. Suggested optimal grafting
density, $\rho_{g}=0.3$, is marked by red.

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 $h^{\dagger}$ ratio, which is found to be
the largest at $\rho_{g}\leqslant 0.3$. One can also see that all the other
fractions, $R_{g}^{\dagger}$, $R_{e}^{\dagger}$, $A^{\dagger}$, and
$B^{\dagger}$, reach their maximum values at $\rho_{g}\approx 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.

Figure 8: (Colour online) Snapshots of the polymer brush structure (shown in
pink) for low and medium densities, $\rho_{g}=0.05$–0.3, at
$T=310\,\textrm{K}$. Note inhomogeneous structure of a brush for
$\rho_{g}=0.05$–0.2.

The narrow interval of grafting densities around $\rho_{g}\approx 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
$T=310\,\textrm{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
$\rho_{g}\approx 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 $\rho_{g}\approx 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 [55, 37] can be invoked. In brief, if the
average separation between the grafted points, $d$, is getting smaller than
the average radius of gyration of an isolated chain, $\bar{R}_{g}$, 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 $d$, assembled normally to the surface,
where the number of monomers $n$ in one blob can be estimated from the
relation $d=l\,n^{\nu}$, where $\nu=0.5882$, $0.5$, and $1/3$ for the good
solvent, $\theta$-point and poor solvent regime, respectively, inside each
blob. The blob interpretation leads to the following scaling law for $\bar{h}$
[37]

$\bar{h}\sim c\,\rho_{g}^{\epsilon}\,N,\qquad\epsilon={\frac{1-\nu}{2\nu}},$
(18)

where $\epsilon\approx 0.35$ for the case of a good solvent within a blob, and
$\epsilon=1/2$ for the case of a $\theta$-point within a blob (desolvated
chain). Similar scaling laws are also expected for both the values of
$\bar{R}_{g}$ and $\bar{R}_{e}$ for each chain in a brush.

Figure 9: (Colour online) Left-hand frame: dependence of the separation $d$
between grafted points of a brush on grafting density $\rho_{g}$. The values
for the radius of gyration of isolated chain, $\bar{R}_{g}$, at
$T=298\,\textrm{K}$ and $310\,\textrm{K}$, based on figure 2, are shown as
horizontal dashed lines. Right-hand frame: the same for the partial
coordination numbers, $n_{p}$ and $n_{w}$, for polymer and water in the
vicinity of each monomer of a chain ($T=298\,\textrm{K}$).

Hence, the changes in the brush height upon increase of $\rho_{g}$ are
governed by two factors. The first one is the excluded volume effect between
the neighbouring chains, that is brought into play when $d(\rho_{g})$
decreases below the $R_{g}$ value. The dependence of $d(\rho_{g})$ on
$\rho_{g}$ is easily established, $d(\rho_{g})=\rho_{g}^{-1/2}$, and it is
shown in the left-hand frame of figure 9. As one can see, at
$T=298\,\textrm{K}$ the value of $d$ drops below $R_{g}$ at $\rho_{g}\approx
0.07$, whereas at $T=310\,\textrm{K}$ this occurs for a much higher grafting
density, $\rho_{g}\approx 0.5$. As the result, at grafting density
$\rho_{g}=0.3$, the brush at $T=298\,\mathrm{K}$ is already deep into the
dense brush regime, whereas at $T=310\,\mathrm{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 $\epsilon$. 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 $n_{p}$ and
$n_{w}$, the partial coordination numbers associated with the polymer-polymer
and polymer-solvent pairs. Their dependencies on the grafting density
$\rho_{g}$ are shown in the right-hand frame of figure 9. Total number of
beads in the first coordination sphere, $n_{p}+n_{w}$, is close to six, as
expected. At a low grafting density, one has strong solvation of polymer
beads, $n_{s}\gg n_{p}$, whereas strong desolvation occurs at high grafting
densities, where $n_{s}\ll n_{p}$. Crossover takes place at the grafting
densities close to $\rho_{g}\approx 0.45$ (see, black dashed line in figure
9). Let us note that at $\rho_{g}=0.3$ (i.e., the supposed optimum grafting
density for the thermoswitchable brush) the polymer chains are still, at least
partially, solvated.

Figure 10: (Colour online) Terms proportional to $\rho_{g}^{\epsilon}$ in the
expressions for $\bar{h}$, $\bar{R}_{g}$, and $\bar{R}_{e}$, see equation
(19), and their fits at fixed values of the exponent $\epsilon=0.35$ and
$0.5$.

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 $N=60=\textrm{const}$ in current simulations, and assume the
presence of some constant terms, $h^{0}$, $R_{g}^{0}$, and $R_{h}^{0}$, in the
scaling expressions

$\bar{h}=h^{0}+C_{h}\,\rho_{g}^{\epsilon},\qquad\bar{R}_{g}=R_{g}^{0}+C_{g}\,\rho_{g}^{\epsilon},\qquad\bar{R}_{e}=R_{e}^{0}+C_{e}\,\rho_{g}^{\epsilon}.$
(19)

Then, we attempted to find the respective intervals for $\rho_{g}$, where the
fits with $\epsilon=0.35$ (solvation regime) and $\epsilon=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:
$h^{0}\approx 2.9$, $C_{h}\approx 19.0$, $R_{g}^{0}\approx 2.2$, $C_{g}\approx
3.6$, $R_{e}^{0}\approx 4.2$, and $C_{e}\approx 13.4$. The intervals being
found are: $0.1<\rho_{g}<0.5$ for the $\epsilon=0.35$ fit and
$0.5<\rho_{g}<1.3$ for the $\epsilon=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 [55, 37, 34, 45].

Therefore, computer simulations of the model thermoresponsive brush considered
in this study reveals that there is some optimum grafting density,
$\rho_{g}\approx 0.3$, at which the ratio between the brush height at
$T=298\,\textrm{K}$ and $310\,\textrm{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 $\rho_{g}=0.3$ and at these
two temperatures. The first one is related to the excluded volume effect
between the neighbouring chains for $\rho_{g}=0.3$. It already promotes the
increase of the brush height at $T=298\,\textrm{K}$
[$d(\rho_{g})<\bar{R}_{g}$], but has no influence on it at $310\,\textrm{K}$
[$d(\rho_{g})>\bar{R}_{g}$], see left-hand frame of figure 9. The second
factor is related to the fact that the monomers are, at least partially,
solvated for $\rho_{g}=0.3$ at $T=298\,\textrm{K}$, see right frame of figure
9. This causes the chains to swell more as compared to the desolvated chains
at $T=310\,\textrm{K}$. With the increase of $\rho_{g}$, the chains at
$T=298\,\textrm{K}$ became desolvated as well, thus decreasing the difference
between the cases of $T=298\,\textrm{K}$ and $310\,\textrm{K}$. With an
increase of the grafting density towards dense brush regime, brush chains
became equally desolvated and stretched both at $T=298\,\textrm{K}$ (below
LCST) and $T=310\,\textrm{K}$ (above LCST), and the difference between these
two temperatures vanishes, see, figure 7.

##  5 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, $\rho_{g}\approx 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 [55, 37]. 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 $\rho_{g}\approx 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 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.

## Acknowledgements

DY and JI thank the CRDF for financial support under Award number
G-202203-68625. DY and OK acknowledge partial support from the NAS of Ukraine
(the grant for research laboratories/groups of young scientists No
07/01-2022(4)). JI thanks to S. Minko, T. Patsahan, and O. Zaichenko for
useful discussions. All authors are grateful to the Armed Forces of Ukraine
for the protection during this research work.

## References

  * [1] Cabane E., Zhang X., Langowska K., Palivan C. G., Meier W., Biointerphases, 2012, 7, No. 1, 9,
doi:10.1007/s13758-011-0009-3.

  * [2] Graziano G., Int. J. Biol. Macromol., 2000, 27, No. 1, 89–97, doi:10.1016/s0141-8130(99)00122-1.
  * [3] Nagase K., Yamato M., Kanazawa H., Okano T., Biomaterials, 2018, 153, 27–48,
doi:10.1016/j.biomaterials.2017.10.026.

  * [4] Cho E. C., Lee J., Cho K., Macromolecules, 2003, 36, No. 26, 9929–9934, doi:10.1021/ma034851d.
  * [5] Ono Y., Shikata T., J. Phys. Chem. B, 2007, 111, No. 7, 1511–1513, doi:10.1021/jp068954k.
  * [6] Maeda Y., Higuchi T., Ikeda I., Langmuir, 2000, 16, No. 19, 7503–7509, doi:10.1021/la0001575.
  * [7] Sun B., Lin Y., Wu P., Siesler H. W., Macromolecules, 2008, 41, No. 4, 1512–1520, doi:10.1021/ma702062h.
  * [8] Ahmed Z., Gooding E. A., Pimenov K. V., Wang L., Asher S. A., J. Phys. Chem. B, 2009, 113, No. 13, 4248–4256, doi:10.1021/jp810685g.
  * [9] Lai H., Wu P., Polymer, 2010, 51, No. 6, 1404–1412, doi:10.1016/j.polymer.2010.01.036.
  * [10] Ono Y., Shikata T., J. Am. Chem. Soc., 2006, 128, No. 31, 10030–10031, doi:10.1021/ja063990i.
  * [11] Philipp M., Kyriakos K., Silvi L., Lohstroh W., Petry W., Krüger J. K., Papadakis C. M., Müller-Buschbaum P., J. Phys. Chem. B, 2014, 118, No. 15, 4253–4260, doi:10.1021/jp501539z.
  * [12] Tavagnacco L., Zaccarelli E., Chiessi E., Phys. Chem. Chem. Phys., 2018, 20, No. 15, 9997–10010,
doi:10.1039/c8cp00537k.

  * [13] Dalgicdir C., van der Vegt N. F. A., J. Phys. Chem. B, 2019, 123, No. 17, 3875–3883,
doi:10.1021/acs.jpcb.9b01644.

  * [14] Tavagnacco L., Chiessi E., Zaccarelli E., Phys. Chem. Chem. Phys., 2021, 23, No. 10, 5984–5991,
doi:10.1039/d0cp06452a.

  * [15] Alaghemandi M., Spohr E., RSC Adv., 2013, 3, No. 11, 3638, doi:10.1039/c3ra22266g.
  * [16] Abbott L. J., Tucker A. K., Stevens M. J., J. Phys. Chem. B, 2015, 119, No. 9, 3837–3845, doi:10.1021/jp511398q.
  * [17] Walter J., Ermatchkov V., Vrabec J., Hasse H., Fluid Phase Equilib., 2010, 296, No. 2, 164–172,
doi:10.1016/j.fluid.2010.03.025.

  * [18] Zanatta M., Tavagnacco L., Buratti E., Chiessi E., Natali F., Bertoldo M., Orecchini A., Zaccarelli E., J. Chem. Phys., 2020, 152, No. 20, 204904, doi:10.1063/5.0007112.
  * [19] Dhamankar S., Webb M. A., J. Polym. Sci., 2021, 59, No. 22, 2613–2643, doi:10.1002/pol.20210555.
  * [20] Abbott L. J., Stevens M. J., J. Chem. Phys., 2015, 143, No. 24, 244901, doi:10.1063/1.4938100.
  * [21] Pérez-Ramírez H. A., Odriozola G., Phys. Chem. Chem. Phys., 2020, 22, No. 32, 17913–17921,
doi:10.1039/d0cp03101a.

  * [22] Shinoda W., DeVane R., Klein M. L., Mol. Simul., 2007, 33, No. 1–2, 27–36, doi:10.1080/08927020601054050.
  * [23] Boţan V., Ustach V. D., Leonhard K., Faller R., J. Phys. Chem. B, 2017, 121, No. 45, 10394–10406,
doi:10.1021/acs.jpcb.7b07818.

  * [24] Bejagam K. K., An Y., Singh S., Deshmukh S. A., J. Phys. Chem. Lett., 2018, 9, No. 22, 6480–6488,
doi:10.1021/acs.jpclett.8b02956.

  * [25] Ilnytskyi J. M., Patsahan T., Sokołowski S., J. Chem. Phys., 2011, 134, No. 20, 204903, doi:10.1063/1.3592562.
  * [26] Kuroki H., Gruzd A., Tokarev I., Patsahan T., Ilnytskyi J., Hinrichs K., Minko S., ACS Appl. Mater. Interfaces, 2019, 11, No. 20, 18268–18275, doi:10.1021/acsami.9b06679.
  * [27] Kim Y., Laradji A. M., Sharma S., Zhang W., Yadavalli N. S., Xie J., Popik V., Minko S., Angew. Chem., 2021, 134, No. 7, e202110990, doi:10.1002/ange.202110990.
  * [28] Soto-Figueroa C., del Rosario Rodríguez-Hidalgo M., Vicente L., Soft Matter, 2012, 8, No. 6, 1871–1877,
doi:10.1039/c1sm07037a.

  * [29] Bautista-Reyes R., Soto-Figueroa C., Vicente L., Modell. Simul. Mater. Sci. Eng., 2016, 24, No. 4, 045004,
doi:10.1088/0965-0393/24/4/045004.

  * [30] del Rosario Rodríguez-Hidalgo M., Soto-Figueroa C., Vicente L., Chem. Phys. Lett., 2018, 695, 170–175,
doi:10.1016/j.cplett.2018.02.016.

  * [31] Lang X., Patrick A. D., Hammouda B., Hore M. J., Polymer, 2018, 145, 137–147,
doi:10.1016/j.polymer.2018.04.068.

  * [32] Groot R. D., Warren P. B., J. Chem. Phys., 1997, 107, No. 11, 4423–4435, doi:10.1063/1.474784.
  * [33] Español P., Warren P., Europhys. Lett., 1995, 30, No. 4, 191–196, doi:10.1209/0295-5075/30/4/001.
  * [34] Descas R., Sommer J. U., Blumen A., Macromol. Theory Simul., 2008, 17, No. 9, 429–453,
doi:10.1002/mats.200800046.

  * [35] Flory P. J., Principles of Polymer Chemistry, Cornell University Press, New York, 1953.
  * [36] des Cloizeaux J., In: Phase Transitions Cargèse, Lévy M., Le Guillou J. C., Zinn-Justin J. (Eds.), 1982, 371–394.
  * [37] de Gennes P. G., Macromolecules, 1980, 13, No. 5, 1069–1075, doi:10.1021/ma60077a009.
  * [38] Ilnytskyi J. M., Holovatch Yu., Condens. Matter Phys., 2007, 10, No. 4, 539, doi:10.5488/cmp.10.4.539.
  * [39] Kalyuzhnyi O., Ilnytskyi J. M., Holovatch Yu., von Ferber C., J. Phys.: Condens. Matter, 2016, 28, No. 50, 505101, doi:10.1088/0953-8984/28/50/505101.
  * [40] Kalyuzhnyi O., Haidukivska K., Blavatska V., Ilnytskyi J., Macromol. Theory Simul., 2019, 28, No. 4, 1900012,
doi:10.1002/mats.201900012.

  * [41] Haydukivska K., Kalyuzhnyi O., Blavatska V., Ilnytskyi J., J. Mol. Liq., 2021, 328, 115456,
doi:10.1016/j.molliq.2021.115456.

  * [42] Haydukivska K., Kalyuzhnyi O., Blavatska V., Ilnytskyi J., Condens. Matter Phys., 2022, 25, No. 2, 23302,
doi:10.5488/cmp.25.23302.

  * [43] Šolc K., Stockmayer W. H., J. Chem. Phys., 1971, 54, No. 6, 2756–2757, doi:10.1063/1.1675241.
  * [44] Guida R., Zinn-Justin J., J. Phys. A: Math. Gen., 1998, 31, No. 40, 8103–8121, doi:10.1088/0305-4470/31/40/006.
  * [45] Binder K., Milchev A., J. Polym. Sci., Part B: Polym. Phys., 2012, 50, No. 22, 1515–1555,
doi:10.1002/polb.23168.

  * [46] Lhuillier D., J. Phys. France, 1988, 49, No. 5, 705–710, doi:10.1051/jphys:01988004905070500.
  * [47] des Cloizeaux J., Phys. Rev. A, 1974, 10, 1665–1669, doi:10.1103/PhysRevA.10.1665.
  * [48] de Gennes P., Scaling Concepts in Polymer Physics, Cornell University Press, London, 1979.
  * [49] Milner S. T., Witten T. A., Cates M. E., Macromolecules, 1988, 21, No. 8, 2610–2619,
doi:10.1021/ma00186a051.

  * [50] Lai P. Y., Binder K., J. Chem. Phys., 1991, 95, No. 12, 9288–9299, doi:10.1063/1.461158.
  * [51] Lai P. Y., Binder K., J. Chem. Phys., 1992, 97, No. 1, 586–595, doi:10.1063/1.463554.
  * [52] Binder K., Eur. Phys. J. E, 2002, 9, No. 3, 293–298, doi:10.1140/epje/i2002-10076-2.
  * [53] Li H., Qian C. J., Sun L. Z., Luo M. B., Polym. J., 2010, 42, No. 5, 383–385, doi:10.1038/pj.2010.9.
  * [54] Grest G. S., Murat M., Macromolecules, 1993, 26, No. 12, 3108–3117, doi:10.1021/ma00064a019.
  * [55] Alexander S., J. Phys., 1977, 38, No. 8, 983–987, doi:10.1051/jphys:01977003808098300.

Ukrainian 0 Моделювання термочутливої полімерної щітки за допомогою
мезоскопічного методу комп’ютерного моделювання Д. Яремчук, О. Калюжний, Я.
Ільницький

Iнститут фiзики конденсованих систем Нацiональної академiї наук України, вул.
Свєнцiцького, 1, 79011, Львiв, Україна