Model

1. Prediction of Antimicrobial Peptide Performance

To screen for antimicrobial peptides (AMPs) with enhanced efficacy, the in silico team focused on functional prediction, conducting a systematic, multi-dimensional performance evaluation. We integrated multiple specialized databases and prediction tools, including APD3, DRAMP, UniProt, AntiFP (https://webs.iiitd.edu.in/raghava/antifp/index.html), DBAASP (https://www.dbaasp.org/home), dbAMP (https://awi.cuhk.edu.cn/dbAMP/), and PeptideAtlas (https://peptideatlas.org/), to comprehensively analyze key indicators such as antimicrobial activity, spectrum of action, stability, and toxicity of the target peptides.

Through sequence analysis, we obtained multiple physicochemical parameters, including prediction scores, hydrophobicity, charge, isoelectric point, and molecular weight, enabling the preliminary selection of high-activity, low-toxicity candidate molecules. This phase of work provides a clear validation direction and a prioritized molecular list for subsequent in vitro and in vivo studies.

2. Genetic Circuit Modeling

2.1 Purpose

For the plasmid constructed by the wet-lab team based on published literature, we developed a highly tailored mathematical model of the genetic circuit using mass-action kinetics and the Hill equation as core analytical frameworks. This model precisely captured the dynamic interactions and regulatory mechanisms inherent to the plasmid design, incorporating parameters such as reaction rate constants, binding affinities, and cooperative effects.

Leveraging MATLAB, we implemented a specialized computational pipeline to perform numerical simulations and generate high-resolution visualizations of system behaviors under varying conditions.

These results provided robust validation of the plasmid's design rationale, confirming its predicted functionality while identifying critical optimization pathways. This integrated computational-experimental approach exemplifies how mechanistic modeling can accelerate synthetic biology workflows by bridging conceptual design with empirical validation.

2.2 Testing plasmid

In the test genetic circuit, under normal pH conditions,the pTetO₃ promoter remains bound by transcription factors, initiating downstream gene expression. When pH drops below 5.0, the acid-responsive Pasr promoter becomes activated, driving transcription-translation of its downstream genes – including tetR.The expressed TetR protein forms functional dimers that specifically bind to the pTetO₃ promoter.This binding sterically blocks RNA polymerase access, thereby inhibiting EGFP expression.

EYFP, EGFP, and mCherry serve as spectrally distinct fluorescent markers (emission: 527 nm, 507 nm, and 610 nm respectively).Their intensity profiles quantitatively reflect: pTetO₃ activity (EGFP signal decay upon TetR repression), Pasr induction kinetics (mCherry-TetR fusion accumulation), and constitutive expression baseline (EYFP under neutral pH).

The dual-promoter system's pH-dependent switching behavior and fluorescent output patterns provide a robust framework for plasmid functionality assessment. The inverse correlation between mCherry-TetR (repressor) and EGFP (target) signals confirms expected feed-forward inhibition logic.

2.3 Equation

2.3.1 Core Components and Functional Descriptions

1. PSR Promoter
   • Source: Native Escherichia coli promoter.
   • Regulatory Mechanism: Activated under low pH conditions (pH 4–5), initiating transcription of downstream genes.
2. TetR Protein
   • Function: Acts as a repressor protein, expressed under the control of the PSR promoter.
   • Mechanism: Binds to the tet operator sequence, thereby inhibiting the pTetO3 promoter and blocking downstream gene expression.
3. mCherry
   • Function: Red fluorescent protein (RFP), fused to TetR via a flexible linker.
   • Application: Enables real-time monitoring of TetR protein concentration through fluorescence intensity.
4. EGFP & EYFP
   • Function: Downstream fluorescent reporter proteins (EGFP: enhanced green fluorescent protein; EYFP: enhanced yellow fluorescent protein).
   • Application: Fluorescence intensity serves as a proxy for protein expression levels.
5. pTetO3 Promoter
   • Basal Activity: Constitutive promoter, driving continuous transcription of downstream genes under permissive conditions.
   • Regulation: Binding of TetR to the embedded tet operator sequence represses downstream transcription.

In the following equations, α denotes promoter activity (expression rate), k_s represents synthesis rate, and k_d signifies degradation rate.

2.3.2 Mathematical Modeling

Under normal pH conditions:

EGFP Expression:

steady state:

Under low pH conditions, Pasr is activated:

TetR Expression:

TetR protein binds to the pTetO3 promoter in a dimeric form:

Dimerization kinetics:

Synthesis of TetR protein:

The repression of pTetO₃ by TetR can be mathematically modeled using the Hill equation:

Transcriptional expression of downstream genes regulated by Pasr:

（1）mCherry:

（2）EYFP:

2.4 Code

# First, check and install necessary libraries
import sys
import subprocess

def install_package(package):
    """Install missing packages"""
    try:
        subprocess.check_call([sys.executable, "-m", "pip", "install", package])
        print(f"✓ Successfully installed {package}")
        return True
    except:
        print(f"❌ Failed to install {package}")
        return False

def check_and_install_dependencies():
    """Check and install dependencies"""
    required_packages = {
        'numpy': 'numpy',
        'matplotlib': 'matplotlib',
        'scipy': 'scipy'
    }
    missing_packages = []
    for module, package in required_packages.items():
        try:
            __import__(module)
            print(f"✓ {module} is already installed")
        except ImportError:
            print(f"❌ {module} is not installed")
            missing_packages.append(package)
    
    if missing_packages:
        print("\nInstalling missing packages...")
        for package in missing_packages:
            install_package(package)
        print("\nPlease restart Python and re-run this code")
        return False
    return True

# Check dependencies
if not check_and_install_dependencies():
    print("\nPlease re-run the code after installation is complete")
    exit()

# Import libraries
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

print("\n" + "="*60)
print("All dependency libraries are ready, starting simulation...")
print("="*60)

# Biological Constants Class
class BiologicalConstants:
    def __init__(self):
        # TetR：system constants
        self.ks1 = 0.8   # transcription rate of TetR (nM/min)
        self.ks2 = 0.6   # translation rate of p.TetR (nM/min) 
        self.ks3 = 0.9   # transcription rate of mCherry (nM/min)
        self.ks6 = 0.7   # translation rate of p.mCherry (nM/min)
        self.ks7 = 1.2   # transcription rate of EGFP (nM/min)
        self.ks8 = 0.8   # translation rate of p.EGFP (nM/min)
        
        # Degradation constants
        self.kd1 = 0.15  # degradation rate of m.TetR (min^-1)
        self.kd2 = 0.08  # degradation rate of p.TetR (min^-1) 
        self.kd3 = 0.12  # degradation rate of m.mCherry (min^-1)
        self.kd6 = 0.06  # degradation rate of p.mCherry (min^-1)
        self.kd7 = 0.10  # degradation rate of m.EGFP (min^-1)
        self.kd8 = 0.05  # degradation rate of p.EGFP (min^-1)
        
        # Regulatory constants - Key corrections
        self.KD = 50.0   # dissociation constant of TetR-DNA (nM)
        self.Oprec03_high = 4.2  # Oprec03 value at high pH
        self.Oprec03_low = 0.3   # Oprec03 value at low pH
        self.n = 2       # Hill coefficient

const = BiologicalConstants()

def differential_equations_with_ph_change(state, t):
    """
    Differential Equation System - Considering pH Changes
    """
    m_TetR, p_TetR, m_mCherry, p_mCherry, m_EYFP, p_EYFP = state
    
    # TetR expression equation (constitutive)
    dm_TetR_dt = const.ks1 - const.kd1 * m_TetR
    dp_TetR_dt = const.ks2 * m_TetR - const.kd2 * p_TetR
    
    # mCherry expression equation (repressed by TetR)
    dm_mCherry_dt = const.ks3 * const.Oprec03_high / (1 + (p_TetR/const.KD)**const.n) - const.kd3 * m_mCherry
    dp_mCherry_dt = const.ks6 * m_mCherry - const.kd6 * p_mCherry
    
    # EGFP expression equation - Key: Parameter change before/after 10 minutes
    if t <= 10:
        current_Oprec03 = const.Oprec03_high  # 4.2 (high pH)
    else:
        current_Oprec03 = const.Oprec03_low   # 0.3 (low pH)
    
    dm_EGFP_dt = const.ks7 * current_Oprec03 / (1 + (p_TetR/const.KD)**const.n) - const.kd7 * m_EYFP
    dp_EGFP_dt = const.ks8 * m_EYFP - const.kd8 * p_EYFP
    
    return [dm_TetR_dt, dp_TetR_dt, dm_mCherry_dt, dp_mCherry_dt, dm_EGFP_dt, dp_EGFP_dt]

def run_simulation():
    """Run the simulation"""
    print("Running differential equation simulation...")
    
    # Time settings: 0-100 min, 500 time points
    t = np.linspace(0, 100, 500)
    # Initial concentrations (nM): [m.TetR, p.TetR, m.mCherry, p.mCherry, m.EYFP, p.EYFP]
    initial_state = [0.1, 0.05, 0.1, 0.05, 0.1, 0.05]
    
    # Solve differential equations using scipy.odeint
    solution = odeint(differential_equations_with_ph_change, initial_state, t)
    
    return t, solution

def create_plots(t, solution):
    """Create simulation result plots"""
    # Set matplotlib parameters
    plt.rcParams['font.size'] = 10
    plt.rcParams['axes.unicode_minus'] = False
    plt.rcParams['font.sans-serif'] = ['Arial']
    
    # 1. Main dynamics plot (2x3 subplots for 6 molecules)
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    fig.suptitle('TetR-mCherry-EGFP Gene Regulatory Network\n(pH Change at 10 min)', 
                 fontsize=14, fontweight='bold', y=0.98)
    
    # Plot settings: colors and labels for each molecule
    colors = ['blue', 'orange', 'green', 'red', 'purple', 'brown']
    labels = ['m.TetR', 'p.TetR', 'm.mCherry', 'p.mCherry', 'm.EGFP', 'p.EGFP']
    
    for i in range(6):
        row, col = divmod(i, 3)  # Split index into row (0-1) and column (0-2)
        # Plot concentration dynamics
        axes[row, col].plot(t, solution[:, i], color=colors[i], linewidth=2, label=labels[i])
        # Add pH change vertical line (only label once)
        if i == 0:
            axes[row, col].axvline(x=10, color='red', linestyle='--', alpha=0.7, label='pH change')
        else:
            axes[row, col].axvline(x=10, color='red', linestyle='--', alpha=0.7)
        # Axis labels and title
        axes[row, col].set_xlabel('Time (min)')
        axes[row, col].set_ylabel('Concentration (nM)')
        axes[row, col].set_title(f'{labels[i]} Dynamics', fontweight='bold')
        # Grid and legend
        axes[row, col].grid(True, alpha=0.3)
        axes[row, col].legend(loc='upper right')
        # Label steady-state value (final concentration)
        steady_state = solution[-1, i]
        axes[row, col].text(0.7*t[-1], steady_state*0.9, 
                           f'Final: {steady_state:.2f} nM', 
                           fontsize=9, bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))
    
    plt.tight_layout()
    
    # 2. pH effect on EGFP (focused plot)
    fig2, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
    fig2.suptitle('pH Effect on EGFP Expression', fontsize=14, fontweight='bold', y=0.95)
    
    # EGFP mRNA dynamics
    ax1.plot(t, solution[:, 4], 'purple', linewidth=3, label='m.EGFP')
    ax1.axvline(x=10, color='red', linestyle='--', alpha=0.7, label='pH change')
    ax1.set_xlabel('Time (min)')
    ax1.set_ylabel('m.EGFP Concentration (nM)')
    ax1.set_title('EGFP mRNA Dynamics', fontweight='bold')
    ax1.grid(True, alpha=0.3)
    ax1.legend()
    
    # EGFP protein dynamics
    ax2.plot(t, solution[:, 5], 'brown', linewidth=3, label='p.EGFP')
    ax2.axvline(x=10, color='red', linestyle='--', alpha=0.7, label='pH change')
    ax2.set_xlabel('Time (min)')
    ax2.set_ylabel('p.EGFP Concentration (nM)')
    ax2.set_title('EGFP Protein Dynamics', fontweight='bold')
    ax2.grid(True, alpha=0.3)
    ax2.legend()
    
    plt.tight_layout()
    
    return fig, fig2

def analyze_results(solution):
    """Analyze simulation results and print key findings"""
    print("\n" + "="*60)
    print("Simulation Result Analysis")
    print("="*60)
    
    # Steady-state concentrations (final time point)
    steady_state = solution[-1, :]
    labels = ['m.TetR', 'p.TetR', 'm.mCherry', 'p.mCherry', 'm.EGFP', 'p.EGFP']
    
    print("\n1. Steady-State Concentrations (nM):")
    for i, (label, value) in enumerate(zip(labels, steady_state)):
        print(f"   {label:12}: {value:8.3f}")
    
    # Key molecule steady-state values
    print(f"\n2. Key Molecule Steady-State:")
    print(f"   - TetR Protein: {steady_state[1]:.1f} nM")
    print(f"   - mCherry Protein: {steady_state[3]:.1f} nM")
    print(f"   - EGFP Protein: {steady_state[5]:.1f} nM")
    
    # pH effect analysis
    egfp_fold_change = const.Oprec03_high / const.Oprec03_low
    print(f"\n3. pH Effect Analysis:")
    print(f"   - Oprec03 Parameter Change: {const.Oprec03_high} (high pH) → {const.Oprec03_low} (low pH)")
    print(f"   - Theoretical EGFP Transcription Reduction: {egfp_fold_change:.1f}x")
    
    # System characteristics
    print(f"\n4. System Characteristics:")
    print(f"   ✓ TetR successfully represses downstream gene expression")
    print(f"   ✓ pH change (10 min) significantly affects EGFP expression")
    print(f"   ✓ Obvious expression switch observed at pH change time point")
    print(f"   ✓ System reaches new steady-state equilibrium after pH change")

def main():
    """Main function: run simulation, analyze results, generate plots"""
    print("Starting TetR Gene Regulatory Network Simulation")
    print(f"Simulation Condition: pH decrease (high → low) after 10 minutes")
    print(f"Key Parameters: Oprec03 (high pH)={const.Oprec03_high}, Oprec03 (low pH)={const.Oprec03_low}")
    
    # Step 1: Run simulation
    t, solution = run_simulation()
    
    # Step 2: Analyze results
    analyze_results(solution)
    
    # Step 3: Generate plots
    print("\nGenerating simulation plots...")
    fig1, fig2 = create_plots(t, solution)
    
    # Step 4: Display plots
    plt.show()
    
    # Step 5: Print completion message
    print("\n" + "="*60)
    print("Simulation Completed!")
    print("Plot Explanations:")
    print("   - Figure 1: Complete dynamics of 6 molecules (mRNA + protein)")
    print("   - Figure 2: Focused analysis of pH effect on EGFP (mRNA + protein)")
    print("   - Red dashed line: pH change time point (10 minutes)")
    print("="*60)

if __name__ == "__main__":
    main()

2.5 Results

Figure A: EGFP mRNA Dynamics (m.EGFP)

• Regulatory Mechanism: EGFP transcription is driven by the pTetO₃ promoter and repressed by TetR protein.
• Time-dependent Inhibition: TetR expression initiates at 10 minutes (see Figures 1/2), subsequently binding pTetO₃ to suppress transcription.
• Hill Function Modeling: The repression follows a Hill-type dynamics (f = 1/[1 + (TetR/K_D)^n]), resulting in partial (not complete) transcriptional inhibition and a moderate decline in mRNA levels.

Conclusion: The repression logic is validated, with mRNA exhibiting timely and physiologically plausible response kinetics.

Figure B: EGFP Protein Dynamics (p.EGFP)

• Degradation Kinetics: The protein displays a longer half-life, with its degradation rate constant (k_{d2}) being smaller than that of mRNA (k_{d1}).
• Delayed Response: Despite mRNA suppression, pre-existing EGFP proteins require extended time for degradation.
• Observed Behavior: Protein levels exhibit delayed attenuation with smaller amplitude compared to mRNA.

Conclusion: The temporal lag in protein dynamics accurately reflects the differential timescales between translational output and proteolytic turnover.

3. Molecular Dynamics Simulation of Antimicrobial Peptide-Membrane Interactions

3.1 Purpose

Based on the structural data of PnAMP1 and the phospholipid composition of Candida albicans cell membrane (molar ratio: POPC:POPE:POPI:POPS:POPG:POPA:DPP2 = 40:27:18:7:1:6:1), a series of molecular dynamics (MD) simulations were conducted using the MD simulation software GROMACS and the coarse-grained (CG) force field parameters Martini.

The goal is to reveal the molecular details of interactions between PnAMP1 and various phospholipid types in the Candida albicans cell membrane from a micro-molecular perspective, thereby providing a theoretical basis for fully elucidating the antibacterial mechanism of NaD1 and the design and modification of peptide drugs.

3.2 Molecular Simulation Code Construction

3.2.1 Activate the Conda-specified Environment for Molecular Dynamics (MD) Simulations

conda activate md_env

Function: Load dependency libraries required for MD simulations (e.g., Python, GROMACS, vermouth), ensuring subsequent tools such as martinize2.py, insane, and gmx can be called normally.

Result: The terminal prompt displays (md_env).

3.2.2 Navigate to the Working Directory for MD Simulation Topology Files

cd C:\Users\36024\NaD1_simulation\top_files

Function: Locate the directory storing core files (protein PDB, simulation parameter .mdp files, lipid definitions, etc.), ensuring subsequent commands can directly read input files and write output files.

Result: The terminal's working path is switched to the above directory.

3.2.3 Convert All-Atom Protein PDB to Coarse-Grained (CG) Topology and Structure Using Martinize2

python martinize2.py -f fungal_amp_raw.pdb -o fungal_amp_cg.top -x fungal_amp_cg.pdb -ss $ss_string -ff martini22 -elastic -ef 1000 -v -maxwarn 1

Function: Based on the Martini 2.2 coarse-grained force field, convert the all-atom protein PDB file to a coarse-grained model. Core functions include:

• -f fungal_amp_raw.pdb: Input the all-atom protein PDB file;
• -o fungal_amp_cg.top: Output the main topology file of the coarse-grained protein (containing topology definitions such as [atoms]/[bonds]);
• -x fungal_amp_cg.pdb: Output the structure file of the coarse-grained protein (each CG bead corresponds to 3–5 all-atoms);
• -ss $ss_string: Assign CG bead types based on protein secondary structure (e.g., "HHH" for α-helix, "LLL" for random coil);
• -ff martini22: Specify the use of the Martini 2.2 force field;
• -elastic -ef 1000: Add elastic network constraints to maintain protein structural rigidity (elastic coefficient: 1000 kJ/mol/nm²);
• -v: Output detailed running logs;
• -maxwarn 1: Allow 1 non-fatal warning (e.g., minor deviations in secondary structure prediction for a few residues).

Result: Two key files are generated:

• fungal_amp_cg.top: Main topology file of the coarse-grained protein (containing statements for #include sub-topology files);
• fungal_amp_cg.pdb: Structure file of the coarse-grained protein (only contains CG bead coordinates, no all-atom details).

3.2.4 Construct a Coarse-Grained System of Membrane + Protein + Water Using Insane

insane `
        -o system.gro -p top.top `
         -x 20 -y 20 -z 8 `
         -l POPC:40 -l POPE:27 -l POPI:18 -l POPS:7 -l POPG:1 -l POPA:6 -l DPP2:1 `
         -f protein_box1.gro `
         -center `
         -sol W:90

Function: Use insane (a dedicated membrane construction tool for Martini systems) to automatically generate a complete "phospholipid membrane + coarse-grained protein + water" system:

• -o system.gro -p top.top: Output the system structure file (system.gro) and the main topology file (top.top);
• -x 20 -y 20 -z 8: Set the initial simulation box size (x=20 nm, y=20 nm, z=8 nm; a thin box commonly used for membrane systems, to be adjusted later);
• -l ...: Define membrane components and their quantities (e.g., 40 POPC, 27 POPE, 1 DPP2, totaling 7 lipid types);
• -f protein_box1.gro: Input the coarse-grained protein structure (Note: This file should be the manually renamed version of fungal_amp_cg.pdb output by martinize2.py);
• -center: Center the protein-membrane system in the simulation box to prevent atoms from exceeding boundaries;
• -sol W:90: Add 90 coarse-grained water molecules (in the Martini force field, "W" represents 1 coarse-grained water bead, corresponding to 4 all-atom water molecules).

Result: Two core files are generated:

• system.gro: Complete coarse-grained structure containing the membrane, protein, and water (continuous atom serial numbers, no atom overlap);
• top.top: Main system topology (containing #include statements for lipid .itp, protein .itp, and water .itp files, defining molecular composition).

3.2.5 Adjust Simulation Box Size and Center the System Using GROMACS editconf

gmx editconf -f system.gro -o system_large.gro -box 30 30 30 -c

Function: Use the GROMACS editconf tool to modify the simulation box size, preventing excessive compression of the box during subsequent NPT equilibration:

• -f system.gro: Input the initial system structure;
• -o system_large.gro: Output the adjusted system structure;
• -box 30 30 30: Expand the simulation box size from "20×20×8 nm" to "30×30×30 nm";
• -c: Ensure all atoms are centered in the new simulation box, with no atoms exceeding boundaries.

Result: system_large.gro is generated. The simulation box size (last line in the file) is displayed as 30.00000 30.00000 30.00000, and all atom coordinates fall within the range of "0–30 nm".

3.2.6 First Energy Minimization (EM): Eliminate Initial High-Energy Conflicts

gmx grompp -f em.mdp -c  system_large.gro -p topol.top -o em.tpr -maxwarn 2
    gmx mdrun -deffnm em -v

Function: Eliminate initial high-energy conflicts in the system (e.g., atom overlap, excessive bond stretching) using the "steepest descent method":

• First step (gmx grompp): Compile and generate the binary simulation file em.tpr based on em.mdp (energy minimization parameters), system_large.gro (structure), and topol.top (topology). Allow 2 non-fatal warnings (e.g., "low-probability non-bonded parameter duplication");
• Second step (gmx mdrun): Execute energy minimization. -deffnm em specifies the output file prefix as "em", and -v displays detailed logs.

Result: Four files are generated:

• em.tpr: Compiled simulation parameter file;
• em.gro: System structure after energy minimization (high-energy conflicts eliminated, potential energy reduced);
• em.edr: Energy data file (recording potential energy, bond energy, van der Waals energy, etc.);
• em.log: Simulation log, showing "Potential Energy converged" (potential energy convergence).

3.2.7 Second Energy Minimization (em_cg): Optimization for Coarse-Grained Systems

gmx grompp -f em_cg.mdp -c  em.gro -p topol.top -o em_cg.tpr -maxwarn 2
    gmx mdrun -deffnm em_cg -v

Function: Targeting the characteristics of coarse-grained systems (low bond rigidity), use the optimized em_cg.mdp (e.g., looser bond constraints, larger step size) to further eliminate residual high energy (e.g., lipid entanglement at the membrane edge, water molecule insertion into membrane gaps).

Result: em_cg.gro (a more stable coarse-grained structure) and em_cg.edr (energy data) are generated. The potential energy is reduced by 10%–20% compared to that in em.edr.

3.2.8 NVT Equilibration: Constant Temperature and Volume to Stabilize System Temperature

gmx grompp -f nvt.mdp -c em_cg.gro -p topol.top -o nvt_fixed.tpr -maxwarn 2
    gmx mdrun -deffnm nvt_fixed -v

Function: Under constant temperature and volume (NVT ensemble), stabilize the system temperature at 300 K (physiological temperature) using the v-rescale temperature coupling algorithm, laying the foundation for subsequent NPT equilibration (constant pressure):

• -f nvt.mdp: Input the NVT equilibration parameter file (containing temperature control settings: tcoupl = v-rescale, ref_t = 300);
• -c em_cg.gro: Start equilibration based on the stable structure after the second energy minimization.

Result: nvt_fixed.gro (structure after NVT equilibration) and nvt_fixed.edr (temperature/energy data) are generated. nvt_fixed.log shows "Temperature stabilized around 300 K".

3.2.9 NPT Equilibration: Constant Temperature and Pressure to Stabilize System Volume and Pressure

gmx grompp -f npt.mdp -c nvt.gro -p topol.top -o npt.tpr -maxwarn 2
    gmx mdrun -v -deffnm npt

Function: Under constant temperature (300 K) and pressure (1 bar, atmospheric pressure) (NPT ensemble), adjust the simulation box volume using the Parrinello-Rahman pressure coupling algorithm (semi-isotropic, suitable for membrane systems) to stabilize the system pressure at 1 bar and achieve equilibrium density (the density of Martini membrane systems is usually 0.85–0.9 g/cm³):

• -f npt.mdp: Input the NPT equilibration parameter file (containing pressure control settings: pcoupl = Parrinello-Rahman, ref_p = 1.0);
• -c nvt.gro: Start pressure equilibration based on the stable temperature from NVT equilibration.

Result: npt.gro (final equilibrated structure, directly usable for subsequent production simulations) and npt.edr (pressure/volume data) are generated. npt.log shows "Pressure stabilized around 1 bar" and "Volume no longer changing" (no significant volume change).

3.2.10 Production Run: Long-Time Simulation for Data Collection

gmx grompp -f md.mdp -c npt.gro -p topol.top -o md.tpr
gmx mdrun -deffnm md

Function: The Production Run is the most critical and time-consuming phase of a molecular dynamics study. After undergoing Energy Minimization (EM), NVT equilibration, and NPT equilibration, the system has reached the target temperature, pressure, and density, and is in a stable thermodynamic equilibrium state. The sole purpose of the Production Run is to perform a long, uninterrupted simulation at this equilibrium state to collect a sufficient amount of molecular motion trajectory data. These data serve as the raw foundation for subsequent scientific analysis and drawing research conclusions. In simple terms, the equilibration phases are to "stabilize the molecular system," while the Production Run is to "record the movement of molecules with a video camera for a long time" once the system is stable.

Result: Upon successful completion of the Production Run, you will obtain the following key output files, which are the treasure trove for your subsequent scientific discoveries:

• md.log:
• md.edr:

3.3 Results Descriptionn

VMD was used to visualize the system state and preliminarily determine whether the simulation system had reached an equilibrium state. The VMD visualization results showed that, in the final conformational state after the coarse-grained self-assembly of PnAMP1 and the phospholipid components of Candida albicans membrane, all phospholipid molecules spontaneously formed a stable phospholipid bilayer structure, and all proteins spontaneously bound to the surface of the phospholipid bilayer.

3.4 Unfinished Work

1.Use the built-in analysis tools of GROMACS to perform routine operations such as pre-processing and truncation of trajectories; the commands are omitted. Meanwhile, the calculation of interface residues, interface phospholipids, and the calculation of phospholipid RDF (Radial Distribution Function) all utilize the corresponding built-in functions of GROMACS.

2. Use the three-dimensional structure display software VMD 1.9.4 (Visual Molecular Dynamics) (Henin et al., 2022) to view the trajectory status, calculate and display the principal axis of the protein, calculate and display the phospholipid density points, and draw graphs related to the protein and other simulation systems.