Phase diagram of two-patch colloids with competing anisotropic and isotropic interactions

TL;DR

This paper maps the phase diagram of two-patch colloids with competing anisotropic and isotropic interactions, revealing various phases and metastable structures relevant to protein aggregation.

Contribution

It provides the first detailed equilibrium phase diagram for two-patch colloids with competing interactions, highlighting phase coexistence and metastability.

Findings

Coexistence of disordered clusters, chains, crystals, and bundles.

Metastability of network bundles versus finite-sized bundles.

Bundles exhibit helical structures similar to sickle cell hemoglobin.

Abstract

Patchy particles are considered to be a good model for protein aggregation. We calculate the equilibrium phase diagram of two-patch colloidal particles undergoing aggregation along with isotropic potential. This investigation demonstrates the coexistence of different phases like disordered clusters, chains, crystals and bundles depending on the relative strength of isotropic and anisotropic potential, these phases are also reported for glucose isomerase proteins. We also show that the formation of network of bundles is metastable against the formation of thermodynamically favored finite sized bundles, which are also metastable along with thermodynamically stable crystals. These bundles appear to be helical in structure similar to that observed in sickle cell hemoglobin.

Equations

$$U({\mathbf{r}}_{i,j},{\mathbf{v}}_i,{\mathbf{v}}_j)=\{\begin{array}{l}\infty r_{i,j}\le \sigma \\ -(u_i+u_a)\sigma <r_{i,j}\le \sigma (1+ϵ)\\ 0r_{i,j}>\sigma (1+ϵ)\end{array}$$

Full text

Human diseases such as cataract, sickle cell anaemia, Alzheimer’s disease, Amyotrophic lateral sclerosis are caused due to undesired protein aggregationsiezen1985opacification ; vekilov2007sickle ; ghiso2002amyloidosis ; patel2015liquid . However, it is also desired to build protein crystals via self assembly to get insight into the structure of proteins from a fundamental point of view berman2006worldwide and for various pharmaceutical applications like targeted drug delivery basu2004protein . Patchy particles are considered to be a very good model for protein aggregation guo2014predictive , they can also be used in imitating atomic orbitals wang2012colloids , tailor made colloidal structures via self assembly manoharan2003dense ; zhang2004self etc. It has already been shown that the anisotropic potential developed by Kern and Frenkel kern2003fluid called patchy particles can mimic certain type of equilibrium protein crystals roberts2014role . Later patchy particles along with isotropic interaction liu2009self ; vekilov2005two ; liu2007vapor ; gliko2005metastable ; gliko2007metastable were found to be closer to experimental results than just isotropic or anisotropic models. It was shown that by playing with the number of patches along with short range isotropic interaction reproduced the liquid-liquid coexistence curves of gamma-crystalline and lyzozyme liu2007vapor proteins. The lyzoyme protein shows a transition from chains to bundles woodard2014gel as seen in amyloid fibers believed to be responsible for some of the neuro-degenerative diseases. Bundle formation is also observed in other biological systems like the actin filaments sept2001thermodynamics . Chain to bundle transition was shown for a two patch particle where irreversible patchy interaction was coupled with reversible isotropic interaction malhotra2018aggregation . Huisman et al.huisman2008phase also observed the formation of bundles and they showed that transition is similar to sublimation transition of polymers. Preisler et al.preisler2013phase also observed the formation of a metastable tubes for the case of particles with a single patch.

In the present work we have simulated two patch colloidal particles interacting via reversible anisotropic and isotropic potential. Tuning both anisotropic and isotropic interaction we are able to observe the existence of a rich phase diagram consisting of phases such as chains + crystals, thermodynamically favored bundles +crystals, disordered clusters + crystals and disordered clusters. Whitelam et al.Whitelam proposed a 2D model where he showed that different pathways for crystallization is possible by tuning the relative strength of the nonspecific and specific interactions. For $2$ patch particles Li et al.li observed the formation of different stable crystal structures when the patch well depth was kept twice of isotropic well depth, but they were not able to observe fiber-bundle transition. Audus et al. Audus also studied the phase behavior for five patch particle with competing isotropic and anisotrpic interaction, where they observed only liquid-liquid phase separation.

Our model consists of hard spheres of diameter $\sigma$, each particle is associated with a unit patch vector $\hat{\bf v}_{i}$ which defines the two oppositely located patches. We have isotropic potential coupled with Kern and Frenkel kern2003fluid anisotropic potential given by

[TABLE]

where $r_{i,j}$ is the distance of centre of mass of the $i_{th}$ particle to the $j_{th}$ particle and $\epsilon$ is the interaction range. $u_{i}$ and $u_{a}$ are the depths of square well for isotropic and anisotropic interactions respectively. $u_{a}$ is non zero only if ${\textbf{\^{r}}}_{i,j}.{\textbf{\^{v}}}_{i}\ >\cos\omega\text{and}\ {\textbf{\^{r}}}_{j,i}.{\textbf{\^{v}}}_{j}>\cos\omega$, where $\omega$ is the semi solid angle of a cone having vertex at the center of a sphere which defines the patch size. In the present work we have used $\sigma=1$, $\epsilon=0.1$ and $\omega=\pi/8$ which ensures that only one bond per patch is possible malhotra2018aggregation . In the present study, all quantities are reported in reduced units, hence bond strength is measured in units of $u/k_{B}T$ where $k_{B}T=1$. Reduced time is given by $t/t_{0}$, where $t_{0}$ is the time taken by a particle to travel it’s own diameter malhotra2018aggregation and the volume fraction is defined as $\phi=\pi/6N_{tot}/L^{3}$. We have kept $\phi=0.02$ which corresponds to a total number of $N_{tot}=4774$ particles. These particles are randomly distributed in a cubic box of fixed length $L=50$ with periodic boundary condition. Every simulation step is divided into two steps, first one is the movement step, where we randomly select $2N_{tot}$ particles which are either rotated with a step size $s_{R}$ or translated with a step size $s_{T}$ prabhu2014brownian . This ensures that every particle undergoes translational and rotational diffusion in an uncorrelated and independent manner in every step. We have fixed $s_{T}=0.013$ and $s_{R}=0.018$ to ensure the correct diffusional behavior prabhu2014brownian . All movement steps which lead to the overlap of spheres or breakage of bonds are rejected. Second step is the cluster construction step, when two monomers are in interaction range and their patches are facing each other then a Patchy $(P)$ bond is formed with probability $\alpha_{a}$ and if a bond already exists it is broken with a probability $\beta_{a}$, such that $\frac{\alpha_{a}}{\alpha_{a}+\beta_{a}}=1-\exp(-u_{a})$ see Fig. 1b. When particles are in interaction range and the patches are not facing each other, they form Non Patchy Isotropic ($NPI$) bond with a probability $\frac{\alpha_{i}}{\alpha_{i}+\beta_{i}}=1-\exp(-u_{i})$ prabhu2014brownian , where $\alpha_{i}$ and $\beta_{i}$ are the probabilities to form and break the isotropic bond respectively see Fig. 1a. The collection of these bonded particles is considered to be a cluster. In the present work we have considered $u_{i}$ ranging from 0 to 2.6 and $u_{a}$ values ranging from 0 to 10.

The phase diagram for two-patch model has been shown in Fig. 2 where we observe the existence of five different phases. At intermediate to high isotropic ($1.39<u_{i}<2.57$) interaction and at low to intermediate anisotropic interaction ($2<u_{a}<6.67$) strength we observe the formation of crystals and $DC$ (disordered clusters) as shown in snapshot of the system in Fig. 3a. As can be observed this region of phase diagram is dominated by isotropic part of the potential and the phase separation is similar to what is observed for system having only isotropic potential babu2006phase . For the well width $\epsilon=0.1$, the isotropic square well system phase separates at $u_{i}\geq 2.57$ babu2006phase while in presence of anisotropic interaction we observe phase separation shifts to $u_{i}<2.57$. Below this region, as we decrease $u_{i}$ and $u_{a}$ we observe the formation of only $DC$ as shown in Fig. 3b. In this region both the isotropic as well as anisotropic part of the potential is relatively small. When $u_{a}$ is increased ($u_{a}>$ 7.69) keeping the same $u_{i}=1.3$ we observe the formation of $DC$ as well as chains see Fig. 3c. This region is dominated by anisotropic part of the potential as $u_{a}>u_{i}$. As mentioned before due to the small patch size, only one bond per patch is possible and as anisotropic bond strength is large, it leads to the formation of chains. Due to the presence of a weak isotropic interaction we also observe a small fraction of $DC$. When we increase $u_{i}$, for the same $u_{a}>7.69$ we observe the formation of crystals along with chains as isotropic part of the potential is playing a competing role (upto $u_{i}=2.14$) as shown in Fig. 3d. As $u_{i}$ is further increased ($u_{i}>2.14$) crystals and bundles along with $DC$ are observed as shown in Fig. 3e. We also observe that the finite bundle that is formed always has helical arrangement. It has already been shown that a bundle with an helical arrangement is a thermodynamically favored state grason2007chirality ; yang2010self ; turner2003twisted , which is probably why we are observing the bundles with crystals Wellems ; Makowski ; Pumphrey ; Jones . For the case of single patch particle it was reported that finite bundles have a lower free energy compared to the lamellar structure which was the equilibrium structure preisler2013phase .

In order to understand the kinetics of formation of different phases, we follow the average number of neighbors bonded via only isotropic interaction $<Z_{NPI}>$ and only patchy interaction $<Z_{P}>$ (in the inset) as a function of reduced time starting from random distribution of particles. For $u_{i}=1.39$ and $u_{a}=3.33$ we observe only disordered clusters where $<Z_{P}>$ and $<Z_{NPI}>$ increases very slowly and equilibrates around a value $10^{-1}$ as the isotropic and anisotropic potential are small. For chains$/DC$ $<Z_{P}>\sim 2$, showing that all the patches are bonded due to high anisotropic interaction strength ($u_{a}=10$) and $<Z_{NPI}>$ equilibrates around $0.17$ due to low isotropic interaction strength $(u_{i}=1.39)$. When $u_{i}=2.3$ and $u_{a}=5$ we observe the formation of crystals as seen by the sudden upturn in the average number of neighbors $<Z_{NPI}>$, which is also observed in the case of $<Z_{P}>$. When the anisotropic part of the potential is increased $(u_{a}=10)$, we also observe that $<Z_{NPI}>$ increases and a sudden upturn is observed similar to crystal structure formation, but $<Z_{P}>\sim 2$ indicates that we are observing chain like structure. A visual inspection of the images as shown in Fig. 4 we observe a network made up of bundles similar to that seen in actin filaments sept2001thermodynamics .

Colloids with short range interaction phase separates into crystal-gas phase and this transition is first order in nature. To know the equilibrium states at different points $(u_{a},u_{i})$ on phase diagram, we have cooled the system from a random configuration to $(u_{a},u_{i})$ as well as heated an already phase separated system (crystallized) configuration to the same point in phase diagram.

In Fig. 5a we have plotted $<Z_{P}>$ as a function of reduced time, the squares indicate the evolution of system to $(3.03,1.39)$ starting from an already phase separated system formed at $(5,2.14)$ (see snapshot to the left of Fig. 5) and triangles indicate the evolution to $(3.03,1.39)$ starting from a random configuration of patchy particles. We observe that the dense crystal phase (indicated by squares) melts completely and meets with the curve starting from a random distribution indicating that the $DC$ formed is the equilibrium phase (see snapshot in the right side of Fig 5a). We also observe that similar trend is observed in the case of $<Z_{NPI}>$ as shown in the inset of Fig 5a. In Fig. 5b squares indicate the evolution of system to $(3.03,2.14)$ starting from a phase separated system formed at $(5,2.14)$ and triangles indicate the evolution of $Z_{P}$ to $(3.03,2.14)$ starting from a random distribution of particles. We observe that dense phase initially evaporates indicated by drop in $<Z_{P}>$ value, but crystal grows again indicating that system formed is metastable in nature (see snapshot to the right top of Fig. 5b). In the inset we can observe that $<Z_{NPI}>$ of the dense phase increases continuously indicating that the crystals are not melting in this particular case. When we start the simulation from a random configuration and quench the system to $(3.03,2.14)$, the system is stuck in a meta-stable state and we observe only $DC$ (see snapshot to the right bottom of Fig. 5b). From a randomly distributed system when we quench the system to $(10,2.14)$ which is in the crystal/bundle region of the phase diagram we observe a network of bundles (see snapshot in Fig. 4). When we start from an already phase separated system formed in the crystal$/DC$ region $(5,2.14)$ we observe that the crystals in the system do not melt, whereas the $DC$ aggregate together to form finite helical bundles as shown in Fig. 3e (also see Supplemental material Supplementary ).

To differentiate between different structures we plot $g(\Theta)$ in Fig. 6, which gives the probability of occurrence of angle $\Theta$ between patch vectors of adjacent particles malhotra2018aggregation . In Fig. 6a we have plotted the $g(\Theta)$ for the crystal structure and for the $DC$ cluster as obtained from the simulations. For the case of disordered cluster we observe that all the angles between $20^{\circ}$ and $140^{\circ}$ are possible and hence these clusters are called as disordered cluster. We have isolated the crystals from the phase separated system and then calculated $g(\Theta)$ for the crystals alone. We observe $2$ prominent peaks for the crystal at $20^{\circ}$ and $160^{\circ}$, which indicates that the constituent spheres of the crystal are aligned along the patch vector. Along with this intermediate angles are also possible with a much smaller probability which may be due to the particles along the surface of the crystal. In Fig.6b we have plotted the $g(\Theta)$ after isolating the bundles and chains from a phase separated system. In this case the anisotropic interaction $u_{a}=10$ is greater than isotropic interaction $u_{i}=1.4$ and we observe chain formation. When the isotropic interaction is higher $u_{i}=2.3$ we observe the formation of bundles. The two peaks we observe are consistent with the fact that for a chain as well as a bundle the patch vector of the spheres is aligned either parallel or anti parallel with the neighboring sphere. It is already shown that once the system forms bundles it will twist into a helical structure and this is considered to be a thermodynamically favored state grason2007chirality ; yang2010self . Similar kind of structures are observed in biological systems like sickle cells turner2003twisted and actin filaments sept2001thermodynamics . Due to the finite twist as predicted by Turner et al. turner2003twisted we are observing only a small tail in the distribution of $g(\Theta)$ at $40^{\circ}$ and $130^{\circ}$ for the parallel and anti parallel patch vectors as shown by arrows in Fig. 6b. Due to the presence of this twist we can easily differentiate between a bundle and a chain apart from the visual confirmation.

In the present work, the formation of various phases such as disordered clusters, chains+crystals and crystals+ thermodynamically favored bundles is achieved by tuning the relative strength of $u_{i}$ and $u_{a}$. The kinetics of formation of helical bundles in our system is distinct from formation of helices as proposed by Chen et al.chen . In their work they studied the helical bundle formation of Janus particles while in the present work we have studied particles with $2$ patch along with isotropic interaction. Chen et al.chen proposed step-by-step addition of individual particles forming small clusters which later fuse together into fibrillar triple helices of finite diameter. While in the present study only when we have relatively high isotropic and anisotropic interaction we observe the presence of thermodynamically favored bundles. We have also shown that by tuning $u_{a}$ and $u_{i}$ we were able to mimic the structures observed for glucose isomerase protein van2018molecular , which include nanorods (chains), crystals, disordered clusters and bundles. The bundles which we observe seems to be helical in structure which is a thermodynamically favored state as reported by Grason et al. grason2007chirality ; yang2010self , which is similar to that observed in sickle hemoglobin fibres turner2003twisted . It will be interesting to study the kinetics as well as the structural properties by varying patch sizes.

We thank HPC padum and Badal of IIT Delhi for providing us the necessary computational resources.