Core-Shell Nanoparticles (examples/core_shell_nps)

This tutorial demonstrates the modeling of a co-assembly system involving a Block Copolymer (PEO-b-PAsp(DET)) and a rigid biological macromolecule (CRISPR-Cas9, PDB: 4OO8).

Key features of this tutorial:

Rigid Body Integration: Importing an external PDB file (4OO8_1_hs.pdb) to represent a complex biological macromolecule (RNP) as a rigid bead in the CG simulation.
Implicit Solvent Modeling: Controlling the assembly behavior (collapse vs. swelling) by tuning the attraction strength scaling factor (\(\alpha\)).
Specific Backmapping: Reconstructing all-atom coordinates while preserving the internal structure of the rigid protein.

Note

Force Field Note: While we assign OPLS-AA parameters here, standard OPLS-AA may require specific tuning for complex proteins/nucleic acids (e.g. OPLS-AA/M). In this context, the force field assignment primarily serves to construct a valid topology for further relaxation or specific force field replacement (e.g., CHARMM, AMBER) later.

Step 1: Chemical Definitions & Rigid Input

We define the block copolymer components (PEO and Asp-derivative) and the rigid protein.

O: PEO unit (Shell-forming block).
N: Functionalized Aspartame unit (Core-forming block).
R: CRISPR-Cas9 RNP, defined by an external PDB file.

import sys
from domd_tools import create_cg_system, build_aa_topology, assign_ff_parameters, get_gmx
import pickle
from misc.logger import logger
logger.setLevel('INFO')

# 1. Define SMILES
smiO = 'OCC' # Polyethylene oxide
smiN = 'OC(=O)C(CC(=O)NCCNCCN)N' # Functionalized Asp monomer (DET)
smiR = 'CC' # Placeholder SMILES for the rigid center

# 2. Define Reaction (Polymerization)
reaction_template = {
    'O-O': {
        'cg_reactant_list': [('O', 'O')],
        'smarts': '[CH3:1].[OH1:2]>>[C:1][O:2]',
        'prod_idx': [0]
    },
    'N-N': {
        'cg_reactant_list': [('N', 'N')],
        'smarts': '[C:1](=[O:2])[OH1:3].[NH2:4][CH1:5]>>[C:1](=[O:2])[N:4][C:5].[O:3]',
        'prod_idx': [0]
    },
    'O-N': {
        'cg_reactant_list': [('O', 'N')],
        'smarts': '[CH3:1].[NH2:2][CH1:3]>>[C:1][N:2][C:3]',
        'prod_idx': [0]
    },
}

# 3. Load Molecules with Rigid Body PDB
mols = {
    'O': {'smiles': smiO, 'file': None},
    'N': {'smiles': smiN, 'file': None},
    # Crucial: 'file' points to the PDB, 'is_rigid' marks it as a rigid body
    'R': {'smiles': smiR, 'file': '4OO8_1_hs.pdb', 'is_rigid': True},
}

Step 2: Generate Initial CG System

We construct a system containing 100 block copolymer chains and 5 rigid RNP complexes.

Block Copolymer: A linear chain of 272 ‘O’ units and 78 ‘N’ units.
Rigid Protein: 5 isolated rigid bodies.

# Define Block Copolymer Topology (O-O-...-O-N-N-...-N)
n_polymer = 100
bonds = [f'{i}-{i+1}' for i in range(272+78-1)]
bondstring = ','.join(bonds)

# Flexible Chain Metadata
meta1 = {
    'type': ['O']*272 + ['N']*78,
    'bond': bondstring,
    'N': n_polymer,
    'is_rigid': False,
    'is_angle': True
}

# Rigid Body Metadata
n_R = 5
meta2 = {
    'type': ['R'],
    'N': n_R,
    'is_rigid': True, # Must match the molecule definition
    'is_angle': True
}

metas = [meta1, meta2]

# Generate System
# Note: iterations=3 helps converge the initial parameter estimation
create_cg_system(mols, reaction_template, metas, boxl=150, iterations=3)

Outputs:

out_chemfast_cg.xml: System topology.
rigid_meta.pkl: Crucial. Contains the mapping between the CG rigid bead and the atoms in the PDB file.

Step 3: Implicit Solvent Assembly (CG Simulation)

This step is performed in the MD engine (HOOMD/GALAMOST). Since we are using implicit solvent, the assembly behavior is controlled by scaling the Lennard-Jones attraction strength (\(\alpha\)).

Tip

\(\alpha\) Tuning of Lennard-Jones Potential The nonbond Lennard-Jones potential is below:

\[U_{LJ}(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \alpha \left( \frac{\sigma}{r} \right)^6 \right]\]

while \(\epsilon\) and \(\sigma\) are determined by the HSP, \(\alpha\) is a user-defined scaling factor to tune the attraction strength. The scaling parameter \(\alpha\) directly correlates with the solvent quality and polymer scaling laws. The empirical guidelines for tuning \(\alpha\) are derived from the study “Brownian dynamics simulation study on the self-assembly of incompatible star-like block copolymers in dilute solution” (Li, B. et al., Phys. Chem. Chem. Phys., 2012, 14, 4964-4970).

\(\alpha < 0.4\) (Hydrophilic/Good Solvent): Weak attraction. The polymer chains swell (\(R_g \sim N^\nu, \nu > 0.5\)).
\(\alpha > 0.6\) (Hydrophobic/Poor Solvent): Strong attraction. The polymer chains collapse (\(R_g \sim N^\nu, \nu < 0.5\)), driving aggregation.

Defining Interactions for Core-Shell Assembly:

For a successful Core-Shell assembly loaded with a rigid cargo (RNP), users are expected to define the interaction potentials based on the specific physicochemical properties of the cargo (e.g., hydrophobicity, charge).

Core Formation: The ‘N’ block (PAsp) should be hydrophobic (\(\alpha > 0.6\)) to form the core.
Shell Stabilization: The ‘O’ block (PEO) should be hydrophilic (\(\alpha < 0.4\)) to stabilize the nanoparticle in solution.
Cargo Loading (Rigid Body ‘R’): * If ‘R’ is hydrophobic or electrostatically attracted to the core, set high affinity between R and N (e.g., \(\alpha_{R-N} \approx 1.0\)). * Ensure ‘R’ has low affinity with the shell ‘O’ (e.g., \(\alpha_{R-O} \approx 0.3\)) to avoid surface adsorption and promote encapsulation.

Protocol: 1. Initialize with out_chemfast_cg.xml. 2. Apply the interaction matrix defined above. 3. Run the simulation until core-shell nanoparticles form and equilibrate. 4. Save the snapshot as CG_pre_eq/final.xml.

Step 4: All-Atom Topology Building (with Rigid Meta)

When backmapping systems with rigid bodies, we must load the rigid_meta.pkl generated in Step 2. This ensures the protein atoms are placed correctly relative to the rigid bead’s center of mass and orientation.

# Path to equilibrated CG snapshot
xmlfile = 'CG_pre_eq/final.xml'

# Load rigid body metadata
rigid_meta = pickle.load(open('rigid_meta.pkl', 'rb'))

# Backmap
# We pass rigid_meta to the builder so it knows how to reconstruct the protein
confs, mol_graphs, box = build_aa_topology(
    mols,
    reaction_template,
    xmlfile,
    chunks_per_d=5,
    rigid_meta=rigid_meta
)

Step 5: Force Field & Output

We use a pure ML strategy here for efficiency on the large polymer systems.

# Assign Parameters using ML strategy for all molecules
strategies = ['ml'] * len(confs)

ffs = assign_ff_parameters(mol_graphs, strategies=strategies)

# Export to GROMACS
get_gmx(mol_graphs, confs, ffs, box)

Final Outputs: The resulting chemfast.gro and chemfast.top will contain the assembled core-shell nanoparticle with the atomistic RNP structure embedded within the PAsp core.