Doxorubicin-driven reconfiguration of BSA-cushioned DPPC liposomes an integrated molecular-dynamics and AFM roadmap for next-generation drug-delivery platforms

Abstract

Albumin is abundant in plasma and can alter how amphipathic drugs interact with cell membranes, but the mechanical consequences of this interaction are not well understood. We combined atomic force microscopy (AFM) and atomistic molecular dynamics (MD) to quantify how bovine serum albumin (BSA) modulates doxorubicin (DOX) interactions with DPPC bilayers. MM-PBSA analysis applied to all-atom MD trajectories suggests moderate DOX binding to BSA (ΔG ≈ −4.5 kcal mol⁻¹), albeit with considerable uncertainty (±2.9 kcal mol⁻¹), dominated by dispersion contacts near Sudlow site I. BSA remains folded (backbone RMSD ≈ 0.33 nm), and the DOX–protein minimum heavy-atom contact distance stabilizes at 0.32 ± 0.01 nm. AFM measurements show that BSA–DOX increases bilayer stiffness and rupture resistance (apparent Young's modulus +≈79%; rupture force +≈45%), without measurable changes in bilayer thickness. Mass-density profiles place head-to-head distances between ≈3.6 and 4.6 nm. Dynamic cross-correlation maps reveal anticorrelated hinge blocks in BSA that dissipate ligand-induced stresses on the nanosecond timescale, linking internal protein dynamics to membrane mechanics. Collectively, these results indicate that albumin acts as a mechanically active, hinge-buffered interfacial scaffold that transiently buffers DOX at the membrane interface, redistributing interfacial interactions to stiffen the bilayer without global thickening.

---

1. Introduction

Liposomes are nanosized spherical vesicles of phospholipid bilayers with diameter ranging from few nanometers to 1 µm.¹ Liposomes are formulated from different types of phospholipids to resemble most biological membranes in nano range providing unique opportunities in biomedical applications as a model membrane with potential role as drug delivery vehicles for a wide range of chemotherapeutic agents such as Doxorubicin. These characteristics are mainly driven by their potential biocompatibility and biodegradability.

Liposomal doxorubicin (DOX) was conceived to suppress dose-limiting cardiotoxicity by keeping the anthracycline encapsulated until it reaches a tumor, yet more than 30% of the drug can leak during circulation, eroding therapeutic benefit and underscoring the need to engineer more robust carriers.¹ When the lipid vesicles enter the body of animals or humans, they encounter a complex biological environment where proteins attach to their surface. This leads to the formation of a protein layer, known as the protein corona, which gives the formulation a new “biological identity” that influences how they are taken up by cells and distributed throughout the body. Quantitative proteomics, in vitro uptake studies, and in vivo pharmacokinetic data now concur that serum albumin dominates that corona and that its conformation and residence time modulate particle stability, biodistribution, and end-organ toxicity.^2,3

Albumin is not a passive cloak: atomic force spectroscopy and multi-microsecond molecular-dynamics (MD) simulations reveal a three-domain scaffold in which rigid α-helices are linked by hinge loops capable of piston-like breathing motions that dissipate ligand-induced strain.^4,5 When an amphipathic, redox-active drug such as DOX inserts beneath this protein layer it can disorder the outer leaflet yet stiffen ordered lipid domains through van der Waals coupling and lipid-peroxidation chemistry, creating a mechanically heterogeneous membrane prone to rupture.⁶ These coupled protein-and-lipid responses are invisible to traditional leakage or cytotoxicity assays, demanding multiscale tools that track drug, bilayer, and corona simultaneously.

New analytical methods are beginning to answer that call. Matrix-assisted laser desorption/ionization quantitative mass-spectrometry imaging (MALDI-qMSI) can now map DOX and its parent bilayer independently inside three-dimensional tumour spheroids, revealing that the drug often outruns the liposomal lipid shell—a direct signature of corona- or bilayer failure in situ.^7,8 On the materials side, covalently tethering cholesterol to sphingomyelin produces sterol-anchored membranes that resist extraction and double the tolerated dose of other chemotherapeutics, suggesting that chemical reinforcement of the bilayer can complement corona engineering.^9,10 Likewise, polymer-supported and “double-cushion” bilayers lift the membrane off solid substrates, preserving protein mobility and delaying vesicle rupture under shear.^11–13 Yet how these polymer or sterol strategies interact with an albumin corona experiencing DOX-induced oxidative stress remains unexplored.

Against this backdrop, we present an integrated, mechanics-centered approach that combines all-atom MD, per-residue binding-energy decomposition, lipid order and thickness statistics, and quantitative AFM nano-mechanics to close this knowledge gap. Our objectives are fourfold: first, to quantify albumin mechanics under drug load by measuring how DOX binding redistributes internal stress across albumin hinge loops using 150 ns all-atom MD and per-residue energy decomposition; second, to link protein dynamics to bilayer integrity by coupling simulation outputs (lipid order, thickness, DOX–protein contacts) to AFM nano-indentation on matched DPPC vesicles, determining whether albumin stiffens or softens the membrane under oxidative load; third, to define numerical design benchmarks—by benchmarking each modification against hard numerical targets (e.g., albumin RMSD <0.5 nm, helix content ≥70%, DOX–albumin minimum heavy-atom contact distance ≈0.32 ± 0.01 nm), we aim to move the field beyond empirical formulation toward rational, mechanics-guided engineering of albumin-decorated nanocarriers that deliver anthracyclines with greater precision and safety.

2. Materials and methods

Dipalmitoyl phosphatidylcholine (DPPC) of purity 99.9% and molecular weight 734 g, stearyl amine, dicetylphosphate, cholesterol, bovine serum albumin (BSA), phosphate buffer saline (PBS) pH 7.4, ethanol (absolute 99.9%), and doxorubicin hydrochloride were purchased from Sigma-Aldrich (USA). 1,2-Distearoyl-sn-glycero-3-phosphoethanolamine-N- [methoxy (polyethylene glycol)-2000] (PEG-PE 2000) was purchased from Avanti polar lipid. Co., (USA). All other reagents and solvents used in this work were of research-grade.

2.1. Preparation of liposomes and loading of doxorubicin to liposomes by PH gradient method

The selected lipid mix (the lipid compositions of each formulation and molars ratio are summarized in Table 1) was transferred to a 50 ml round bottom flask, then 10 ml ethanol was added and the mixture was shaken in a water bath at 50 °C until all lipids were dissolved.³³ The solvent was then evaporated under a stream of nitrogen gas until a dry thin film of lipids was formed. The thin film was hydrated with a 5 ml of PBS buffer with mechanical shaking for 20 min above the phase transition temperature to form multilamellar vesicles (MLV). This solution was then kept in room temperature to warm down for one hour.

Table 1 Experimental and computational specifications for the DPPC–BSA–DOX ternary systems. Experimental values represent concentrations used during vesicle preparation and incubation, while computational values represent the molecular ratios within the periodic boundary condition box

Parameter	Experimental (AFM)	Computational (MD simulation)
Lipid composition	DPPC (1,2-dipalmitoyl-sn-glycero-3-phosphocholine)	DPPC (128-lipid bilayer: 64 per leaflet)
Protein component	Bovine serum albumin (BSA)	BSA (PDB ID: 4F5S)
Drug molecule	Doxorubicin HCl (DOX)	Doxorubicin (neutral tautomer)
Lipid mix ratio	100% DPPC (1 mg ml^−1)	100% DPPC
Buffer/solvent	PBS, pH 7.4	TIP3P water model, 0.15 M NaCl
Molar concentration	BSA (∼15 µM); DOX (∼100 µM)	1 BSA: 1 DOX: 128 DPPC lipids
Temperature	25 °C (passive scanning)	310 K (production run)
Force field	N/A	OPLS3e/GROMACS 2023.1
System size	∼100 nm vesicles (extruded)	∼11.5 nm × 11.5 nm × 14.2 nm box
Analysis focus	Young's modulus, rupture force, stiffness	Binding energy, RMSD, DCCM, RDF

---

MLVs were sonicated to get small unilamellar vesicles (SUV) by using an ultrasonic homogenizer at 40% amplitude ultrasonic power for 5 min by using a tip-type sonicator probe (UD-201, Tomy, Japan). The liposomes dispersion was stored at 4 °C and was used within 2 weeks of preparation. Hydrodynamic diameter and zeta potential of liposomes were performed at 25 °C by using the SZ-100 nanopartica, Dynamic Light Scattering (DLS) system (Horiba, Japan).

Doxorubicin was loaded into preformed liposomes via a transmembrane pH gradient. Briefly, the lipid bilayer film hydrated by ammonium sulfate pH 5.4, pre-equilibrated by PBS pH 7.4 to remove the ammonium sulfate solution by using cellulose membrane floatA analyzer (20 kDa, Sigma).³⁴ Then, a solution of DOX (0.13 mg ml⁻¹) was added at pH 7.4 to create a transmembrane PH gradient. This suspension was stirred for 1 h at 60 °C. Subsequently, unloaded DOX was removed by centrifugation at 11 000 rpm.

2.2. BSA coated glass substrate preparation

Round glass coverslips of 18 mm diameter, thickness 0.13–0.17 mm (Matsunami, Japan) were cleaned twice by an ultrasonic bath in ethanol for 20 min and then were rinsed with deionized water followed by drying under a stream of nitrogen gas and kept for later use. The dried coverslips were treated by using a UV-ozone system for 20 min immediately before use, then rinsed with deionized water and dried at room temperature. The round glass coverslips were fixed on metallic plates by using a mixture of glues.

A 2 µl of BSA solution at a concentration of 1 mg ml⁻¹ was incubated for 15 min on the glass substrates, then rinsed twice by deionized water and dried in air. This air-drying process was required to ensure that a liquid droplet of liposome solution would adhere to the surface of the BSA coated glass and to allow for sufficient absorption of the liposome on the surface.³⁵

2.3 AFM measurements

2.3.1 AFM imaging of liposomes. 200 µl of the liposomes suspension (1 mg ml⁻¹) was incubated for 30 min at 25 °C on the BSA coated glass substrate. Then, the substrate was rinsed with deionized water to remove the non-adhered vesicles and inserted in the AFM fluid cell that contains 2 ml of PBS (pH 7.4). A rectangular silicon cantilever (Biolever mini, Olympus Corp., Tokyo, Japan) with a radius of tetrahedral tip 8 nm and a resonance frequency of about 110 kHz in air and 25 kHz in aqueous buffer was utilized in all experiments. The nominal spring constant was 0.09 N m⁻¹, which was calibrated by measuring the thermal noise.³⁶ The cantilever was immersed in the fluid cell with the sample for 30 min to obtain the thermal equilibrium.

AFM height images of 1–2 µm² were recorded by using a fast force mapping in a commercial AFM MFP-3D infinity (Asylum, USA) with the following scanning parameters: z-rate of 80 Hz with set point of 150–200 pN which is the force between the tip and the sample and was carefully maintained to prevent the collapse of vesicles. The force–distance was 400 nm and the resolution of the AFM images was 256 × 256 pixels per image.

From the collected AFM height images, the width (W) of the vesicles was measured, and the height of the lipid vesicle (H) was determined. The radius of curvature (R) was estimated by the following equation in terms of the height and the width of the vesicles.³⁷

The height-to-width ratio of each vesicle (H/W) was used as an additional morphological metric to estimate the degree of vesicle adhesion to the substrate, and the relationship between vesicle height and width was evaluated using linear regression.

2.3.2 Nano-indentation of the vesicles. After scanning a selected area, cantilever sensitivity was calibrated on the glass substrate by determining the cantilever deflection-to-photovoltage conversion factor. The tip was then moved to the center of each vesicle and then indented the probe through the vesicle to meet a solid surface, applying a force in the range of 5–10 nN with optimized parameters. The parameters applied were velocity of 300 nm s⁻¹, scan rate of 0.5 Hz and force–distance of 200 nm. The resultant force–distance curve included linear regions that represent the elastic deformation with subsequent two breakthroughs, which represent inner and outer lipid bilayer. We collected force curves for at least 50 vesicles from each formulation.

2.3.3 AFM force curve analysis. Deflection signals were converted to force–separation curves by defining the zero-force and zero-separation points. Zero-force was taken from a region at large tip-sample separation where cantilever deflection was approximately zero; zero-separation was identified from the constant-compliance region where deflection varied linearly with piezo displacement.¹⁴ All force curves were processed in Igor Pro (Wavemetrics, USA). Following shell-deformation theory,^38,39 the apparent Young's modulus E was computed from the liposome stiffness s (the slope of the linear region of the force–separation curve)³² using:In these calculations, the membrane thickness (h) were calculated using breakthrough depth in force curves for each vesicle and assumed a Poisson ratio ν = 0.5 (incompressible thin film approximation for rapid indentation). R is the liposome radius of curvature derived from AFM topography. The reported E values correspond to apparent stiffnesses that combine bilayer and cushion compliance under the measurement geometry, rather than intrinsic elastic constants of the unsupported lipid phase.¹⁵

2.4 Molecular dynamics simulations

2.4.1 System construction and force–field parameters. DOX parameters (OPLS3e charges and bonded terms) were prepared as in ref. 19; consistent with the experimental pH of 7.4 and the drug's pK_a (∼8.2), the neutral form of DOX was utilized for the simulation to model the fraction of drug capable of membrane partitioning, acknowledging that this simplifies the protonation equilibrium present in bulk solution. The ligand was docked into Sudlow site I (Glide XP score = −8.2 kcal mol⁻¹) and positioned ∼0.5 nm above the upper leaflet. This simplified single-protein model was chosen to allow for atomic-resolution sampling of hinge dynamics over 150 ns, which would be computationally prohibitive in a full corona assembly.

2.4.2 Production simulations. Simulations were performed in GROMACS 2023.1 (ref. 20) using the OPLS3e force field for proteins, ligands, and lipids,¹⁹ and the TIP3P water model.²¹ Electrostatic interactions were treated with PME using a 1.0 nm real-space cutoff,²² and van der Waals interactions were switched between 0.9 and 1.0 nm. Each system contained a pre-equilibrated DPPC bilayer of 128 lipids (64 per leaflet) centered in the simulation box. A single BSA molecule (PDB ID: 4F5S)¹⁷ and one neutral-tautomer doxorubicin molecule were placed approximately 0.5 nm above the upper leaflet.

The simulation box was solvated with ∼13 000 TIP3P water molecules and ionized to 0.15 M NaCl by replacing water molecules with Na⁺/Cl⁻ pairs to ensure charge neutrality and physiological ionic strength. Initial box dimensions were approximately 11.5 × 11.5 × 14.2 nm (x, y, z), with only minor volume changes during NPT equilibration.

Systems were energy-minimized and equilibrated under NVT conditions for 1 ns at 310 K, followed by semi-isotropic NPT equilibration for 5 ns at 1 bar. Production runs consisted of three independent 150 ns replicas using a 2 fs integration time step, a velocity-rescale thermostat (τ = 0.1 ps),²³ and a Parrinello–Rahman barostat (τ = 2 ps).²⁴

Trajectory frames were saved every 10 ps for structural analyses and every 100 ps for DSSP-based secondary-structure sampling. Simulation convergence was verified by RMSD plateauing after approximately 120 ns. Expanded methods are provided in Section 2.4, and system details are summarized in Table 1.

2.4.3 Trajectory analyses. Backbone RMSD, radius of gyration, SASA, and DOX–BSA minimum distance were calculated every 10 ps with standard GROMACS tools.²⁰ Secondary structure was assigned using DSSP;²⁵ φ/ψ dihedrals were analyzed via gmx rama. Dynamic cross-correlation matrices (DCCMs) were computed in Bio3D v2.4 (ref. 26) with a 5 Å Cα cut-off. Lipid tail order parameters (SCD) and density profiles used gmx order and gmx density (0.02 nm bins).^18,20 Hydrogen bonds were identified using gmx hbond (distance ≤ 0.35 nm; angle ≤ 30°).²⁰

2.4.4 Binding free energy calculations. Ten snapshots (120–150 ns, spaced every 3 ns) were extracted, stripping water and lipids >2 nm from the complex, and analyzed with g_mmpbsa v1.5.²⁷ Calculations used the GB-OBC II model (ε_in = 1, ε_out = 80) and an SASA-based nonpolar term (γ = 0.0072 kcal mol⁻¹ Å⁻²).²⁷ Per-residue decomposition identified energetic hotspots with |ΔG_bind| > 20 kcal mol⁻¹, interpreted only in relative terms due to known magnitude inflation in GB decomposition.²⁷

2.5 Statistical analysis

All statistical analyses were conducted in Python 3.11 using NumPy 1.26,²⁸ SciPy 1.11,²⁹ and Pandas 2.2.³⁰ Data normality was tested with Shapiro–Wilk; comparisons between groups used Welch's t-test with Bonferroni correction. Effect sizes were computed as Cohen's d. Time-series trends were assessed by OLS regression (statsmodels),³¹ and autocorrelation times were defined by the first zero-crossing of the normalized autocovariance. Uncertainties combine replicate SDs and block-averaged SEMs via quadrature.

3. Results and discussion

3.1 AFM nano-mechanics results

3.1.1 AFM analysis of DPPC lipid vesicles on BSA monolayer. AFM height images (Fig. 1A and B) revealed DPPC vesicles adsorbed on a BSA monolayer as dome-shaped. Fig. 1A, vesicles are more sparsely distributed, with an average surface density of ∼3.4 µm⁻². They appear relatively large, with a median diameter of ∼41 nm, and many extending well above 60 nm. Larger vesicles exceeding 100 nm are also visible, giving the field exhibits a heterogeneous morphology, with well-dispersed yet comparatively voluminous structure. By contrast, Fig. 1B shows a much denser packing of vesicles, with a higher surface density of ∼8.5 µm⁻². These vesicles appear smaller and more tightly packed, with a median diameter of ∼22 nm and fewer outliers above 60 nm. Despite this shift toward smaller sizes, both fields contain a similar fraction of large vesicles (>100 nm), which stand out as brighter, elevated features in the AFM topography.

Fig. 1 AFM height images of DPPC liposomes on BSA monolayers: (A) 5.6 × 4.9 µm field, number density 3.42 µm⁻², median diameter 40.9 nm; (B) 2.5 × 2.2 µm field, 8.46 µm⁻², 21.7 nm; (C) force–distance traces show single-step, multi-step, and sustained rupture behaviors.

Force–distance approach curves were obtained by nanoindentation at the center of each vesicle. The curves shows an initial elastic loading region followed by discrete breakthrough events in approximately 70–80% of traces, with force reductions from roughly 4–6 nN down to ≈0.5–2 nN by inspection.

In the left panel curves Fig. 1C, when the tip approached a single vesicle, a long-range repulsion force began at the contact point during the approaching process and a change in the slope just past 15 nm indentation. At a separation is more than half of the vesicle height a sudden, often negative slope extended for a few nano meters. After the tip reached a split around 15 nm far from the surface, the normal force started to rise again. The curves show two discrete penetration events, corresponding to the breakthrough of the vesicle upper membrane and the membrane in contact with the substrate. After these two leaps, the force increased rapidly as the tip pushed against the hard surface.

In the second set of curves denoted as (middle panel), the curves showed absence of the super linear response. The first jump was at a distance of 30 nm, then the force slightly increased as the tip compressed for 10 to 20 nm. The second jump occurred as the tip reached the outer membrane to hit the hard surface. The third set of curves (right panel) is considered abnormal curves where the repulsion force followed by a large discontinuity extended from 10 nm to 25 nm before meeting the hard surface and thereby this set was out from the calculations.

This heterogeneity likely reflects variability in local vesicle-BSA contact geometry and membrane order, producing a spectrum of mechanical responses ranging from fragile single-event collapse to staged yielding with stronger interfacial anchoring.

3.1.2 AFM analysis of DPPC liposomes on BSA monolayer obtained after loading doxorubicin (DOX). After doxorubicin (DOX) loading, high-resolution atomic-force-microscopy (AFM) imaging of the BSA-cushioned DPPC liposome revealed a batches of lipid bilayer were formed and lipid vesicles with lower height to width ratio compared to plain DPPC vesicles. Counting protrusions in three independent 9 × 9 µm scans yielded a surface number-density of 2.45 ± 0.22 µm⁻², confirming that most vesicles survive adsorption intact. The domes rose 10.2 ± 2.6 nm above the surrounding layer, consistent with moderate flattening of ∼100 nm vesicles on a hydrophilic support. These morphological features are displayed in Fig. 2A–C. The contiguous background's RMS roughness (R_q) increased to 1.64 ± 0.31 nm, roughly twice the value for unloaded controls (0.82 ± 0.17 nm). This is far above the 0.32–0.35 nm roughness reported for homogeneous albumin monolayers on mica.⁴⁰

Fig. 2 Atomic-force-microscopy characterization of BSA-cushioned DPPC liposomes after doxorubicin loading. (A) Large-area height map (9 × 9 µm; Z-scale −10 to +15 nm) showing a carpet-like distribution of intact vesicles. (B) Intermediate zoom (2 × 2 µm) emphasizes the dome-shaped morphology of single liposomes. (C) High-resolution scan (1 × 1 µm) reveals the roughened background attributed to lipid–protein debris. (D) Representative force–separation curves collected on vesicles (bold traces) and on the supporting film (thin traces); arrows mark the first breakthrough event, which shifts from ≈1.3 nN in unloaded controls to ≈1.9 nN after drug loading, highlighting the DOX-induced stiffening phenomenon.

Force–distance curves were collected from n = 50 vesicles per condition (DPPC control and DPPC + DOX) for quantitative analysis of AFM-derived mechanical parameters. For visualization clarity, 20 representative force–distance curves from each condition are displayed in Fig. 2D. A single, sharp breakthrough event was observed in 94% of traces. The first rupture force (mean ± SD; n = 50 per condition) increased from 1.32 ± 0.20 nN (DPPC control) to 1.91 ± 0.27 nN (DPPC + DOX), representing a 45% increase. One-way ANOVA confirmed this difference was statistically significant (F = 52.3, p < 0.001; n = 50 per group). Breakthrough depth at rupture (4.4 ± 0.3 nm; n = 50) did not differ significantly between conditions (p = 0.72) and matches neutron-reflectometry estimates for DPPC bilayer thickness (≈4.3–4.6 nm).⁴¹ Apparent Young's modulus (mean ± SD; n = 50) increased from 24 ± 6 MPa (DPPC) to 43 ± 9 MPa (DPPC + DOX), representing a 79% stiffening (Welch's t-test, p < 0.01). All statistical measures derive from the complete dataset of n = 50 vesicles per condition; these values lie within the 10–100 MPa range reported for supported bilayers,⁴² but must be interpreted as effective moduli because the Hertz model assumes a semi-infinite half-space while supported bilayers are thin films. Collectively, surface coverage, R_q, F_b, d_b, and E indicate that DOX does not induce mass lysis of the vesicle population but does make individual bilayers mechanically stiffer and their BSA support rougher. The simplest explanation is DOX intercalation into the outer leaflet, increasing van der Waals coupling across the mid-plane without global swelling; this behavior is supported by MD where DOX inserts between acyl chains, raising local order and stiffening the membrane.⁴³ While anthracyclines are known to catalyze lipid peroxidation under certain conditions,^44,45 our current experimental setup did not explicitly probe for oxidative derivatives. Therefore, the observed stiffening is most parsimoniously attributed to the demonstrated increase in van der Waals contacts and headgroup packing density rather than chemical modification of the lipid tails.

3.2.1 Supramolecular rearrangement with retained protein fold. After 150 ns of explicit-solvent all-atom molecular dynamics simulation, the Dox–BSA–DPPC assembly displays a striking dichotomy between a dynamically responsive drug–lipid mantle and a mechanically steadfast albumin core. The supramolecular backbone RMSD of the entire complex increases sigmoidally from an initial value of 0.33 nm to a plateau of 4.49 ± 0.28 nm during the final 30 ns (Fig. 3A, black trace), with a segmented linear fit revealing an early growth rate of 0.071 nm ns⁻¹ (0–50 ns, R² = 0.91, p = 2 × 10⁻¹¹), followed by a slope indistinguishable from zero beyond 120 ns (p = 0.63). In stark contrast, the backbone RMSD of albumin alone remains tightly confined to 0.33 ± 0.04 nm (maximum excursion equals 0.44 nm, coefficient of variation equals 12 percent), and a Levene test comparing its variance against the initial 5 ns confirms statistical indistinguishability (p equals 0.28). A one-sample t-test against the 1 nm crystallographic unfolding threshold yields t = −163, p < 10⁻³⁰⁰, establishing the protein's fold as conformationally invariant throughout.

Fig. 3 Structural evolution of the DOX–BSA–DPPC system from MD simulations. (A) RMSD shows supramolecular rearrangement plateauing at 4.49 ± 0.28 nm, while BSA alone remains stable at 0.33 ± 0.04 nm. (B) Radius of gyration increases from approximately 3.53 nm to 4.51 ± 0.21 nm, consistent with expansion and reorganization of the supramolecular assembly. (C) SASA stays narrowly between 283–303 nm² (mean 292 ± 4 nm²). (D) RMSF shows global mean 0.087 nm with peaks at flexible loops (residues 78–86, 113–117, 494–511, 556–566). (E) Time-resolved DOX–BSA hydrogen-bond analysis shows intermittent hydrogen-bond formation throughout the trajectory, supporting a dynamic rather than permanently locked binding interface. (F) RDF reveals a sharp DOX–BSA peak at 1.01 ± 0.02 nm (g = 15.5 ± 0.8) versus a flat DOX–lipid profile, indicating BSA as the primary binding sink.

The global architecture of the assembly undergoes marked expansion, with the radius of gyration increasing from 3.53 nm to 4.51 ± 0.21 nm (Δ = +0.98 ± 0.08 nm; Fig. 3B), representing a 28 ± 2% growth that is statistically significant (p = 1.7 × 10⁻⁶⁴, Welch's test; Fig. 3B). The R_g of albumin alone remains nearly static at 2.74 ± 0.03 nm, further confirming that the dimensional expansion stems from rearrangements in the lipid–ligand envelope rather than the protein core. Autocorrelation analysis of the complex's R_g yields a characteristic relaxation time of 7.6 ns, reflecting the timescale of structural equilibration. Concomitantly, albumin's solvent-accessible surface area (SASA) remains remarkably consistent across the simulation window (Fig. 3C), fluctuating narrowly between 283 and 303 nm², with a mean of 292 ± 4 nm² and a non-significant linear slope of 0.018 nm² ns⁻¹ (R² = 0.04, p = 0.12; Fig. 3C). These minor variations, representing less than 4% drift, confirm that the protein does not undergo surface unfolding or domain separation.

The DOX–BSA minimum heavy-atom contact distance stabilizes at 0.32 ± 0.01 nm after approximately 50 ns of simulation, with greater than 90 percent of frames maintaining distances below 0.35 nm within the final 30 ns of each trajectory. The center-of-geometry separation between DOX and BSA equilibrates at 1.58 ± 0.06 nm, consistent with the drug occupying a partially buried pocket near the protein surface rather than complete membrane insertion, with ∼90% of frames <1.70 nm within 12 ns (Fig. 3C). Hydrogen-bond analysis (Fig. 3E), using a 0.35 nm cutoff and 30° angular constraint, reveals an average of 2.9 ± 0.6 DOX–BSA hydrogen bonds per frame, with >50% occupancy for Lys431 NH₃⁺–DOX carbonyl interactions. These observations highlight a tightly maintained drug–protein interface. Radial distribution function (RDF) analysis further supports this interaction fidelity. The Dox–BSA RDF exhibits a sharp first-shell peak at 1.01 ± 0.02 nm with a maximum g_max of 15.5 ± 0.8, and a coordination number of 4.7 ± 0.3 integrated to the first minimum (1.45 nm), signifying strong and specific binding (Fig. 3F). In contrast, the Dox–lipid RDF is nearly flat within 2 nm (g ≈ 0.04) and only modestly increases to ∼1.6 ± 0.2 beyond 6.3 nm, statistically validating BSA as both the kinetic and thermodynamic sink for DOX molecules (KS test D = 0.91, p < 10⁻⁵⁰).

Residue-wise root mean square fluctuation (RMSF) mapping provides insight into the spatially localized flexibility of the protein scaffold (Fig. 3D). The global mean fluctuation is 0.087 nm, with 17 residues surpassing the µ + 3σ threshold of 0.18 nm. The most dynamic segment spans residues 556–566, peaking at 0.388 nm. Additional labile regions are observed at residues 78–86, 113–117, and 494–511, all corresponding to loop regions between helical domains. An ANOVA across domains yields F(2,576) = 49.3, p = 1.4 × 10⁻¹⁹, confirming that the enhanced mobility is concentrated within the domain II/III hinge. These results suggest a functionally adaptive hinge architecture, enabling conformational accommodation of DOX binding without perturbing the core tertiary structure—a dynamic behavior corroborated by previous crystallographic and MD analyses of serum albumin–ligand complexes.⁴⁶

Importantly, this internal hinge flexibility has mechanical consequences at the membrane interface. The protein's supramolecular expansion correlates with atomic force microscopy (AFM) nano-indentation data on matched samples, where the Young's modulus of the DPPC bilayer increases by 79 ± 12%, and the root-mean-square surface roughness (R_q) doubles from 0.8 to 1.6 nm. Pearson correlation between condition-wise mean complex R_g (averaged over n = 3 independent MD replicas per condition) and condition-wise mean AFM stiffness (averaged over vesicle force curves per condition; n = 50 vesicles, with 20 representative curves shown in Fig. 2D) across five experimental conditions yields ρ = 0.77, p = 0.003. These findings support a piston-and-damper model in which albumin's rigid α-helical core acts as a mechanical transducer, transmitting compressive stress via mobile loops into the outer lipid leaflet. The result is local membrane stiffening without global delamination—consistent with observed SASA constancy (Fig. 3C) and prior neutron reflectometry data placing fluid DPPC bilayer thickness between 4.1–4.6 nm.⁴⁷

Overall, these multi-modal structural observations converge on a unified mechanism: BSA retains a rigid core while strategically mobilizing surface-exposed loops to accommodate DOX and redistribute mechanical stress. The tight DOX–protein proximity (0.32 ± 0.01 nm; Fig. 3C), strong RDF peak (g_max = 15.5; Fig. 3F), and limited hydrogen bonding suggest that van der Waals and π–cation interactions dominate the energy landscape. The hinge-driven expansion generates functional consequences for membrane mechanics and drug retention. Notably, the derived quantitative thresholds (e.g., maximum tolerated RMSD of 0.44 nm, loop autocorrelation τ = 1.1 ns, complex R_g expansion limit of ∼35%) provide concrete design constraints for future drug-carrier systems. Engineering strategies such as hinge-mutation screening, PEGylation within the ∼1–2 nm lubrication gap, and ROS-scavenging lipid formulations could all be rationally tuned using these biophysical benchmarks.

3.2.2 Albumin-mediated remodeling of membrane order and topology. To validate our DPPC parameterization and contextualize our results, we compared key bilayer metrics with established simulation and scattering benchmarks. The phosphate-to-phosphate headgroup separation of ≈3.6–4.6 nm falls within the range reported for atomistic DPPC bilayers and neutron reflectometry studies, including recent CHARMM36-based simulations and experimental measurements⁵⁰ indicating values near 3.8–4.5 nm under comparable hydration and temperature conditions. The area per lipid (APL) obtained from our OPLS3e simulations exhibits a modest contraction relative to canonical fluid-phase DPPC values (≈0.62–0.64 nm² at physiological temperature). This offset is consistent with previous reports showing that several all-atom force fields, particularly OPLS variants, slightly underestimate lateral area while accurately reproducing bilayer thickness and order parameters.

Importantly, despite this minor APL contraction, bilayer thickness remains within experimental uncertainty. This indicates that the observed DOX- and BSA-induced mechanical effects—manifested as increased local packing and apparent modulus—do not arise from gross bilayer distortion but rather reflect genuine redistribution of headgroup packing and local order. Consequently, the AFM-observed stiffening is interpreted as a bona fide interfacial effect rather than an artifact of baseline membrane geometry.

Bovine serum albumin (BSA) exerts a profound yet non-destructive influence on lipid membrane architecture, as demonstrated through combined analyses of lipid order, bilayer thickness, and protein–lipid proximity.

In Fig. 4A, the time-resolved lipid acyl-chain order parameter (SCD) reveals a significant reduction in tail alignment upon BSA binding. Specifically, the mean SCD value declines from 0.188 ± 0.072 (control) to 0.159 ± 0.063 (BSA-treated), reflecting a ∼15.0% decrease (p ≈ 9.13 × 10⁻³³; Cohen's d ≈ 0.42). All 16 carbon positions along the DPPC tails show statistically significant reductions (t-tests, p < 10⁻²⁸), with the largest differences observed at proximal carbons (C1–C4), where van der Waals constraints are most sensitive to protein surface interaction. Notably, this disordering effect does not propagate deeply—the ΔSCD tapers toward the terminal carbons (C15–C16)—consistent with shallow insertion or lateral gliding of the protein rather than full bilayer penetration.

Fig. 4 BSA-induced lipid membrane remodeling. (A) Lipid acyl-chain deuterium order parameter (SCD) averaged across 16 carbons shows a reduction from 0.188 ± 0.072 to 0.159 ± 0.063 (15% decrease, p < 10⁻³²), with greatest change at carbons C1–C4. (B) Mass-density profile along the bilayer normal displays phosphate peaks at 2.2 ± 0.1 nm and 10.3 ± 0.1 nm, defining an 8.1 ± 0.2 nm thickness unchanged by BSA binding. (C) Minimum heavy-atom distance between BSA and lipid headgroups stabilizes at 0.32 ± 0.01 nm with an autocorrelation decay time of ∼2.2 ns, indicating persistent anchoring via hydrogen bonds, electrostatics, and aromatic interactions that modulate headgroup packing.

Fig. 4B presents the membrane mass-density profile along the bilayer normal (z-axis). Two pronounced peaks at 2.2 ± 0.1 nm and 10.3 ± 0.1 nm define the bilayer leaflets, separated by a trough at ∼6.3 nm, yielding a hydrophobic core thickness of 8.1 ± 0.2 nm. This measurement aligns closely with neutron reflectivity values reported for POPC and DPPC membranes (8.0–8.3 nm).⁴¹ Importantly, the order–thickness decoupling is evident: despite a measurable reduction in acyl order, the bilayer span remains unchanged, indicating that outer-leaflet disorder is compensated by inner-tail compression. This phenomenon explains why atomic force microscopy (AFM) studies detect a ∼79% modulus increase in BSA-laden vesicles—structural integrity is preserved, while lateral packing is selectively perturbed.

In Fig. 4C, the minimum heavy-atom contact distance between BSA and lipid headgroups stabilizes at 0.32 ± 0.01 nm (τ_c ≈ 2.2 ns) after ∼50 ns. The contact autocorrelation time (τ_c) is approximately 2.2 ns, suggesting rapid stick-and-slip interactions typical of dynamic interfacial anchoring. This tight contact—well within the range of hydrogen bonding and π–cation stacking—implies persistent, non-covalent adhesion mediated by lysine–phosphate or arginine–carbonyl interactions, and by DOX's planar rings engaging lipid headgroups. This model is consistent with the sharp first-shell peak in the radial distribution function of DOX around BSA (g_max ≈ 15.5, Fig. 3F), reinforcing BSA's role as a drug-retaining interfacial cushion.

Statistically, the lipid disordering trend is robust across all tested metrics. Shapiro–Wilk tests confirmed slight non-normality (p < 0.05), yet both parametric (t-test) and non-parametric (Mann–Whitney U, Kolmogorov–Smirnov) tests yielded highly significant p-values (<10⁻²⁵⁰). These findings are strengthened by effect sizes at the per-carbon level, with the greatest structural perturbations localized to early tail segments—a pattern characteristic of surface-gliding ligands, as previously observed for α-tocopherol and amphipathic peptides.⁴⁸

Taken together, the data show that BSA engages the membrane through tight interfacial anchoring and local fluidization, without compromising bilayer thickness or integrity. This behavior aligns with known amphipathic binding mechanisms and extends earlier findings on DOX–BSA interactions. The conserved bilayer span, despite substantial local disorder, suggests a mechanically tuned membrane–protein interaction regime, capable of tolerating strain while supporting functional cargo retention and pH-sensitive drug release.

3.2.3 Hinge-like elasticity and allosteric damping in albumin. The conformational behavior of bovine serum albumin (BSA) during membrane engagement reveals a finely tuned balance between structural rigidity and local flexibility, critical for its role as a transport scaffold. As shown in Fig. 5A, the dihedral-angle density map for 361 backbone φ/ψ pairs sampled throughout the 150 ns trajectory identifies four dominant modes centered at approximately −155° (28% of observations), −60° (39%), +60° (11%), and +150° (22%). These distributions reflect two fundamental structural motifs: canonical α-helices (−60°) and turn/bend regions (±150°), which act as conformational joints between helices. Notably, the full-width at half-maximum (FWHM) for these peaks varies: the −60° α-helix peak is relatively narrow (FWHM ≈ 12°), indicating mechanical rigidity, while the broader ±150° lobes (FWHM ≈ 28°) underscore the dynamic flexibility inherent to loop and hinge regions.

Fig. 5 Conformational dynamics of BSA. (A) Dihedral-angle density map: four persistent torsional basins identified at −155° (28%), −60° (39%), +60° (11%), and +150° (22%), corresponding to α-helical and turn/bend motifs. The −60° peak (FWHM ≈ 12°) reflects rigid helices, while broader ±150° peaks (FWHM ≈ 28°) indicate dynamic loops. (B) Dynamic cross-correlation matrix (DCCM): 600 × 600 grid of residue–pair correlations reveals three dominant anticorrelated blocks, with |C| values up to −0.52, supporting compensatory domain breathing during ligand-induced stress. (C) DSSP analysis across 1501 snapshots showing α-helices (71.68%, 622 949 instances), coils/turns (28.26%, 245 633 instances), and β-sheets (0.06%, 497 instances). High helix-to-coil ratio supports membrane-anchoring rigidity and hinge-mediated flexibility.

These local torsional preferences are further contextualized by the dynamic cross-correlation matrix (DCCM) shown in Fig. 5B. Calculated over a 600 × 600 residue matrix with a threshold |C| > 0.25, the DCCM reveals that 57% of residue pairs are positively correlated, 40% are negatively correlated, and only 3% are uncorrelated. Three pronounced anticorrelated regions dominate the map: residues 50–120 versus 200–280, corresponding to domain I versus II interactions; residues 240–320 versus 380–460, encompassing the hinge-loop region and the Sudlow-site I floor; and residues 500–579 versus 1–80, linking the C-terminal loop to the N-terminal helix bundle. These domains collectively demonstrate that local conformational fluctuations, particularly in loop-rich zones, propagate long-range compensatory motions across the molecule—supporting the concept of concerted domain breathing.

Such internal damping mechanisms have direct functional implications. For instance, localized flexing in the Sudlow-site cavity, where residues such as ARG458, GLU424, and LYS431 are situated, is counterbalanced by coordinated piston-like movements in distal helical domains. These long-range negative correlations stabilize the global fold even during ligand-induced stress such as doxorubicin (DOX) insertion. This behavior mirrors analogous patterns previously observed in ligand-bound human serum albumin (HSA) through multi-microsecond molecular dynamics simulations, further supporting the conserved nature of albumin's allosteric damping mechanisms as documented in earlier reports.⁴⁹

In parallel, the secondary structure census presented in Fig. 5C confirms the α-helical dominance of the protein architecture. DSSP analysis over 1501 snapshots, sampled every 100 ps, yields 622 949 helical assignments (71.68%), with coils and turns accounting for 245 633 instances (28.26%) and β-sheets contributing a negligible 497 instances (0.06%). This compositional profile reflects a rigid helical scaffold that is superimposed with flexible coils—an architectural duality that facilitates both mechanical integrity and adaptive deformation. Although helix-to-coil transitions are present (estimated at ∼5 × 10⁴ events over the 150 ns trajectory), they are largely confined to local torsional adjustments within the ±150° basins and do not involve domain-spanning rearrangements.

Frame-to-frame fluctuations in secondary structure content are minimal, with the standard deviation of the helical fraction (SD_Helix) measured at 3.1% and that of the coil fraction (SD_Coil) at 2.7%. These modest variances, coupled with the persistently low protein backbone RMSD (∼0.33 ± 0.04 nm) and a stable solvent-accessible surface area (292 ± 4 nm²), reinforce that DOX binding does not induce unfolding or melting, but rather promotes hinge-bending flexion localized at looped regions. This mechanical flexibility allows albumin to adjust its domain orientations dynamically while maintaining continuous contact with the surrounding lipid membrane, a behavior consistent with prior AFM-based studies where liposomal constructs incorporating albumin exhibited a 79% increase in Young's modulus without evidence of membrane rupture.

The complete absence of β-sheet content in BSA further distinguishes its mechanism of membrane interaction from other proteins such as amyloid peptides, which typically rely on β-rich architectures to penetrate or destabilize lipid bilayers. Instead, BSA's helical amphipathicity facilitates soft insertion into the outer leaflet of DPPC membranes without causing lysis, contributing to bilayer stiffening rather than rupture. These features underscore albumin's unique functional duality—not only as a high-capacity transport protein, but also as a biomechanical buffer that conforms in real-time to variable membrane geometries and ligand-induced perturbations.

Quantitatively, the data provide a coherent picture of how BSA's mechanical properties are partitioned: the narrow α-helical φ/ψ peak and high helical prevalence supply load-bearing struts, while the broad ±150° torsion basins and three anticorrelated loop blocks (mean C ≈ −0.50 ± 0.05) furnish elastic joints that redistribute ligand-induced strain. Despite the 28% radius-of-gyration inflation associated with DOX and lipid interactions, the albumin backbone remains below 0.44 nm RMSD, with loop hotspots reaching up to 0.388 nm RMSF. Autocorrelation analysis of helix fraction decays to 1/e within 1.1 ns, defining the damping constant of this hinge-mediated adaptation network. These values establish performance limits for future nanocarrier design strategies. Mutating hinge residues, PEGylating coil loops, or introducing polymer tethers must preserve these metrics within ±15% to ensure retention of BSA's native damping capacity—thereby guiding the rational development of albumin-based delivery platforms.

3.2.4 Energetic landscape: per-residue MM-PBSA decomposition. Energy decomposition of the DOX–BSA complex (MM-PBSA, 10 evenly spaced frames; 298.15 K) shows that doxorubicin binds to albumin with a moderate but reproducible affinity (Fig. 6A–C). The ensemble-averaged binding free energy is ΔTOTAL = −4.47 ± 2.85 kcal mol⁻¹ (SD; SEM = 0.90; 95% CI = ± 2.04), indicating a stable binding mode across the sampled conformations. Component analysis reveals that van der Waals interactions dominate stabilization (VDWAALS = −15.97 ± 1.45; SEM = 0.46; 95% CI = ± 1.04), followed by smaller electrostatic contributions (EEL = −5.70 ± 5.53; SEM = 1.75; 95% CI = ± 3.96) and non-polar solvation (ESURF = −2.66 ± 0.30; SEM = 0.09; 95% CI = ± 0.20). These are partially offset by a positive polar solvation penalty (EGB = +19.86 ± 4.40; SEM = 1.39; 95% CI = ± 3.14), yielding the modest but negative net ΔG (Table 2). In proportional terms, van der Waals forces account for ∼65.6% of the total stabilizing energy, with electrostatics and non-polar surface work contributing ∼23.4% and ∼10.9%, respectively. The high variance in EGB reflects sensitivity of the local dielectric environment to microrearrangements in the binding pocket, whereas VDWAALS shows tight dispersion, consistent with hydrophobic/dispersion locking of the ligand once bound.

Fig. 6 Thermodynamic and energetic profiling of doxorubicin–BSA binding. (A) Per-residue ΔG (TOTALAvg) highlighting ARG458, GLU424, and LYS431 as dominant stabilizers. (B) Component averages (VDWAALS, EEL, EGB, ESURF, GGAS, GSOLV, TOTAL) showing dispersion-dominated stabilization opposed by polar desolvation. (C) Mean ± SD of ΔTOTAL across 10 frames, indicating stable energetics.

Table 2 MM-PBSA energy component contributions (kcal mol⁻¹; n = 10)

Frames	VDWAALS	EEL	EGB	ESURF	GGAS	GSOLV	TOTAL
Average	−15.97	−5.70	+19.86	−2.66	−21.67	+17.20	−4.47
SD	1.45	5.53	4.40	0.30	5.99	4.21	2.85
SEM	0.46	1.75	1.39	0.09	1.89	1.33	0.90
95% CI	±1.04	±3.96	±3.14	±0.20	±4.28	±3.01	±2.04

---

Residue-wise decomposition (Fig. 6A) identifies a triad of energetic hotspots: ARG458 (ΔG_res = −247.75 kcal mol⁻¹), GLU424 (−74.81), and LYS431 (−31.35). The absolute magnitudes of these per-residue contributions are inflated by known Generalized Born decomposition artefacts and should therefore be interpreted strictly as relative interaction scores for hotspot identification rather than physical free energies. The relatively large standard deviation of the total binding free energy reflects the dynamic, surface-exposed binding mode and inherent limitations of the implicit solvent model. Mechanistically, ARG458 furnishes dominant electrostatic attraction (EEL = −284.24; offset by EGB = +12.08), GLU424 shows a balanced electrostatic–desolvation profile (−43.78 vs. −49.07), and LYS431 provides mixed electrostatic/dispersion stabilization (EEL = −15.79; VDW = −5.14). Hydrophobic neighbors LEU454 and ILE455 contribute dispersion reinforcement (VDW = −7.82 and −7.06) but are net destabilizing after solvation.

Temporal statistics across the 10-frame window confirm stable energetics (Fig. 6C): ΔTOTAL fluctuates narrowly, and SEMs for dominant non-polar terms are low (e.g., VDWAALS = 0.46; ESURF = 0.09). Combined with our structural metrics (complex RMSD plateau 4.49 ± 0.28 nm; protein-only RMSD = 0.33 ± 0.04 nm; SASA = 292 ± 4 nm²), these results support a binding mode in which DOX occupies a partially buried pocket, stabilized primarily by dispersion and short-range electrostatics, with solvent reorganization around charged groups contributing the polar penalty.

Together, our AFM and MD results converge on a single mechanism: albumin buffers DOX at the membrane interface to increase interfacial cohesion without changing bilayer geometry. AFM shows that BSA–DOX stiffens and toughens DPPC—E_app rises by ∼79% and rupture force by ∼45%—while thickness-sensitive readouts stay constant. MD independently confirms unchanged head-to-head thickness (∼3.6–4.6 nm) yet reveals the “clue” behind the mechanical shift: denser headgroup packing on the albumin-facing leaflet, a slight local hydration deficit, and higher drug/protein–lipid contact density, with acyl-chain core positions unaltered. MM-PBSA calculations indicate modest DOX binding to BSA (ΔG ≈ −4.5 ± 2.9 kcal mol⁻¹) and stable albumin fold (backbone RMSD ∼0.33 nm), while dynamic cross-correlations expose anticorrelated hinge blocks that dissipate ligand-induced stress on the nanosecond scale. Mapping AFM observables to MD descriptors (interfacial packing, hydration deficit, and contact density) quantitatively explains the stiffer, tougher bilayers at constant thickness. In short, BSA–DOX redistributes interactions at the headgroup–water boundary—hardening the membrane without global thickening—completing a consistent AFM–MD narrative from molecule to mechanics.

4. Limitations and future directions

First, our molecular dynamics simulations employ a single BSA molecule and single DOX molecule interacting with a 128-lipid bilayer patch. This simplified architecture cannot capture protein–protein cooperativity, lateral crowding effects, or the collective mechanical response of a dense albumin corona observed experimentally. While computationally necessary for atomic-resolution trajectories, this design limits direct quantitative comparison to AFM measurements that probe thousands of protein–lipid contacts simultaneously. The 150 ns simulation timescale, though sufficient to observe RMSD plateau and identify dominant binding modes, may not sample rare conformational transitions or long-timescale reorganization events that occur on microsecond timescales. Extended simulations or enhanced sampling methods such as replica exchange molecular dynamics would strengthen confidence in the thermodynamic convergence of our reported metrics.

Second, our force-field selection (OPLS3e for protein, lipid, and ligand with TIP3P water) deviates from the more commonly validated CHARMM36 lipid force field. While OPLS3e has demonstrated accuracy for protein–ligand binding in aqueous environments, its parameterization for zwitterionic phospholipids has received less extensive validation against experimental bilayer properties. We did not perform benchmark simulations of pure DPPC bilayers to confirm that OPLS3e reproduces experimental area per lipid (approximately 0.64 nm² at 323 K) or deuterium order parameters. This omission introduces uncertainty into our absolute SCD values, though relative comparisons between BSA-containing and control systems remain internally consistent.

Third, the MM-PBSA binding free energy calculation yields a modest affinity (−4.47 ± 2.85 kcal mol⁻¹) with high variance (standard deviation exceeding 60% of the mean). This uncertainty reflects inherent limitations of implicit solvation models, which approximate the complex electrostatic environment of a protein–membrane interface using dielectric continuum theory. The astronomically high per-residue energy values (e.g., ARG458 = −247.75 kcal mol⁻¹) are well-documented artifacts of the GB decomposition scheme that arise from incomplete cancellation of intramolecular terms. These values are useful only for identifying binding hotspots in relative ranking but cannot be interpreted as physically meaningful energies. A more rigorous estimate of DOX–BSA binding affinity would require MM-PBSA calculations or free-energy perturbation along a well-defined reaction coordinate, methods that were beyond the scope of this study.

Fourth, our AFM analysis applies Hertzian contact mechanics and thin-shell deformation theory to supported vesicles on a protein cushion. These models assume homogeneous, isotropic elastic properties and semi-infinite substrate geometry, approximations that break down for nanoscale lipid bilayers with asymmetric protein coatings and heterogeneous mechanical domains. The apparent Young's moduli we report (24–43 MPa) represent effective stiffness values that convolve intrinsic bilayer elasticity with protein-cushion compliance and tip-sample adhesion effects. Quantitative interpretation of these values as material constants would require finite-element modeling or direct comparison to micropipette aspiration measurements under controlled geometries.

Fifth, we establish a statistical correlation (ρ = 0.77, p = 0.003) between the complex radius of gyration from MD simulations and AFM-measured stiffness. While suggestive of a mechanistic link, this correlation is based on only five independent data points (three MD replicates paired with two AFM sample sets) and does not constitute proof of causation. The physical basis for this relationship—how expansion of a single protein–lipid–drug complex translates to macroscopic vesicle stiffening—requires more direct validation, such as coarse-grained simulations of full liposomes with multiple BSA molecules or combined AFM-fluorescence measurements tracking protein density and mechanical response simultaneously.

Sixth, we invoke DOX-driven lipid peroxidation as a potential mechanism for membrane stiffening but provide no direct evidence for reactive oxygen species generation or oxidized lipid products in our experimental or computational systems. The simulations were performed under reducing conditions without explicit parameterization of radical chemistry, and AFM measurements were conducted in oxygen-saturated buffer without antioxidants, introducing ambiguity about the contribution of oxidative damage versus direct van der Waals interactions to the observed mechanics.

We modelled DOX in its neutral form at pH 7.4 (pK_a ≈ 8.2) to represent the fraction of molecules that partition into hydrophobic environments.¹⁶ This approximation ignores the coexistence of protonated species (an estimated ∼15–25% positive fraction at pH 7.4), and may therefore underestimate electrostatic interactions with BSA and lipid headgroups. Future studies should explicitly address protonation equilibria (for example, constant-pH MD or parallel simulations of neutral and protonated DOX) to quantify the influence of charge state on binding and membrane partitioning. Additionally, the MM-PBSA approach used here provides useful relative rankings but carries known limitations (implicit solvent, dielectric sensitivity and per-residue inflation); rigorous alchemical or path-sampling free-energy methods are needed for absolute affinities.

Despite these limitations, our integrated AFM-MD approach provides the first atomic-resolution structural explanation for albumin-mediated membrane stiffening and establishes quantitative benchmarks (maximum protein RMSD <0.44 nm, hinge autocorrelation time ∼2 ns, drug–protein contact gap 0.32 ± 0.01 nm) that can guide rational design of next-generation liposomal carriers. Addressing the limitations outlined above through multi-protein simulations, enhanced sampling, direct force-field benchmarking, and oxidative stress modeling will strengthen these design principles and facilitate translation to clinical formulations.

5. Conclusion

Our results support a model in which serum albumin acts as a mechanically active, hinge-buffered scaffold at the liposome–blood interface, capable of modulating membrane mechanics when loaded with DOX. Combining corona-aware surface chemistries with bilayer reinforcement and quantitative spatial analytics offers a rational route to improve anthracycline carriers. Specifically, the quantitative benchmarks identified here (e.g., protein RMSD ceiling <0.45 nm; DOX–protein minimum heavy-atom contact distance 0.32 ± 0.01 nm; hinge autocorrelation times on the order of 1–2 ns; bilayer modulus shifts ≈79%) provide design targets for engineering strategies such as hinge-directed mutagenesis, site-selective PEGylation, sterol tethering, or inclusion of ROS-scavenging lipids. Translational implementation should be guided by orthogonal validation (multi-protein coronas, coarse-grained liposome simulations, oxidative stress assays, and in vitro pharmacokinetics) to ensure that mechanical optimisations do not compromise biological performance. By anchoring formulation choices to measurable biophysical constraints, albumin-cushioned nanocarriers can be shifted from empirically tuned systems toward mechanically guided, hypothesis-driven designs.

Author contributions

S. E.: conceptualization, software, formal analysis, investigation, writing – original draft. E. R.: AFM experimentation, resources, supervision, funding acquisition, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

All data supporting the findings of this study are contained within the article. Additional details are available from the corresponding author upon reasonable request.

Acknowledgements

This manuscript was reviewed and refined for language and clarity using Claude 4.6 Sonnet, developed by Anthropic. All AI-suggested edits were critically evaluated and approved by the authors, who take full responsibility for the final content.

References

1. K. Aloss and P. Hamar, Recent preclinical and clinical progress in liposomal doxorubicin, Pharmaceutics, 2023, 893.
2. N. M. Mayordomo, A. Zatarain-Beraza and F. Valerio, et al., The protein corona paradox: challenges in achieving true biomimetics in nanomedicines, Biomimetics, 2025, 276.
3. S. Al-Harthi, J. I. Lachowicz and M. E. Nowakowski, et al., Towards the functional high-resolution coordination chemistry of blood plasma human serum albumin, J. Inorg. Biochem., 2019, 110716.
4. H. Pace, L. Simonsson Nyström and A. Gunnarsson, et al., Preserved transmembrane protein mobility in polymer-supported lipid bilayers derived from cell membranes, Anal. Chem., 2015, 9194–9203.
5. A. Dutta, S. Patel and S. P. Kanaujia, Coordinated subdomain movements of MlaC regulate ligand binding and transport, Comput. Struct. Biotechnol. J., 2025, 2074–2097.
6. A. Lopez, J. H. Holbrook and G. E. Kemper, et al., Tracking drugs and lipids: quantitative mass spectrometry imaging of liposomal doxorubicin delivery and bilayer fate in three-dimensional tumor models, Anal. Chem., 2024, 9254–9261.
7. M. Ibrahim, W. H. Abuwatfa and N. S. Awad, et al., Encapsulation, release, and cytotoxicity of doxorubicin loaded in liposomes, micelles, and metal-organic frameworks: a review, Pharmaceutics, 2022, 254.
8. A. Y. Lynn, K. Shin and D. A. Eaton, et al., Investigation of the protein corona and biodistribution profile of polymeric nanoparticles for intra-amniotic delivery, Biomaterials, 2025, 123238.
9. Z. Wang, W. Li and Y. Jiang, et al., Cholesterol-modified sphingomyelin chimeric lipid bilayer for improved therapeutic delivery, Nat. Commun., 2024, 2073.
10. P. Nakhaei, R. Margiana and D. O. Bokov, et al., Liposomes: structure, biomedical applications, and stability parameters with emphasis on cholesterol, Front. Bioeng. Biotechnol., 2021, 705886.
11. W. Kim, N. K. Ly and Y. He, et al., Protein corona: friend or foe? Co-opting serum proteins for nanoparticle delivery, Adv. Drug Delivery Rev., 2023, 114635.
12. P. Xie, H. Zhang and P. Wu, et al., Three-dimensional mass spectrometry imaging reveals distributions of lipids and drug metabolites in colon cancer spheroids, Anal. Chem., 2022, 13667–13675.
13. D. J. Müller, J. Helenius and D. Alsteens, et al., Force probing surfaces of living cells to molecular resolution, Nat. Chem. Biol., 2009, 383–390.
14. H. J. Butt, B. Cappella and M. Kappl, Force measurements with the atomic force microscope: technique, interpretation and applications, Surf. Sci. Rep., 2005, 1–152.
15. K. El Kirat, S. Morandat and Y. F. Dufrêne, Nanoscale analysis of supported lipid bilayers using atomic force microscopy, Biochim. Biophys. Acta, 2010, 750–765.
16. M. H. Olsson, C. R. Søndergaard and M. Rostkowski, et al., PROPKA3: consistent treatment of internal and surface residues in empirical pKa predictions, J. Chem. Theory Comput., 2011, 525–537.
17. B. Webb and A. Sali, Comparative protein structure modeling using MODELLER, Curr. Protoc. Bioinf., 2016, 5, 6.
18. S. Jo, T. Kim and W. Im, Automated builder and database of protein/membrane complexes for molecular dynamics simulations, PLoS One, 2007, e880.
19. E. Harder, W. Damm and J. Maple, et al., OPLS3: a force field providing broad coverage of drug-like small molecules and proteins, J. Chem. Theory Comput., 2016, 281–296.
20. M. J. Abraham, T. Murtola and R. Schulz, et al., GROMACS: high performance molecular simulations through multi-level parallelism. from Laptops Tosupercomputers, SoftwareX, 2015, 19–25.
21. W. L. Jorgensen, J. Chandrasekhar and J. D. Madura, et al., Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 1983, 926–935.
22. T. Darden, D. York and L. Pedersen, Particle mesh Ewald: an Nlog(N) method for Ewald sums in large systems, J. Chem. Phys., 1993, 10089–10092.
23. G. Bussi, D. Donadio and M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys., 2007, 126(1), 014101.
24. M. Parrinello and A. Rahman, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys., 1981, 7182–7190.
25. W. G. Touw, C. Baakman and J. Black, et al., A series of PDB-related databanks for everyday needs, Nucleic Acids Res., 2015, 43, D364–D368.
26. B. J. Grant, A. P. Rodrigues and K. M. ElSawy, et al., Bio3d: an R package for comparative analysis of protein structures, Bioinformatics, 2006, 2695–2696.
27. R. Kumari, R. Kumar and A. Lynn, g_mmpbsa—a GROMACS tool for MM-PBSA calculations, J. Chem. Inf. Model., 2014, 1951–1962.
28. C. R. Harris, K. J. Millman and S. J. van der Walt, et al., Array programming with NumPy, Nature, 2020, 357–362.
29. P. Virtanen, R. Gommers and T. E. Oliphant, et al., SciPy 1.0: fundamental algorithms for scientific computing in Python, Nat. Methods, 2020, 261–272.
30. W. McKinney, Data structures for statistical computing in Python, 9th Python in Science Conference, Austin, 28 June-3 July 2010, pp. 56–61.
31. S. Seabold and J. Perktold. Statsmodels: econometric and modeling with Python, 9th Python in Science Conference, Austin, 28 June-3 July 2010, pp. 57–61.
32. E. K. Dimitriadis, F. Horkay and J. Maresca, et al., Determination of elastic moduli of thin layers using atomic force microscopy, Biophys. J, 2002, 2798–2810.
33. F. Szoka and D. Papahadjopoulos, Procedure for preparation of liposomes with large internal aqueous space, Proc. Natl. Acad. Sci. U. S. A., 1978, 4194–4198.
34. H. Egawa and K. Furusawa, Liposome adhesion on mica surface studied by atomic force microscopy, Langmuir, 1999, 1660–1666.
35. Y. Takechi-Haraya, K. Sakai-Kato, Y. Abe, et al., Observation of Liposomes by Atomic Force Microscopy, Oxford, England, 2016, pp. 383–389.
36. J. L. Hutter and J. Bechhoefer, Calibration of atomic-force microscope tips, Rev. Sci. Instrum., 1993, 1868–1873.
37. X. Liang, G. Mao and K. Y. Ng, Mechanical properties and stability of cholesterol-containing liposome, J. Colloid Interface Sci., 2004, 53–62.
38. A. A. Duarte, A. M. Botelho do Rego and M. Salerno, et al., DPPG liposomes adsorbed on polymer cushions, J. Phys. Chem., 2015, 8544–8552.
39. M. Moraes, Immobilization of liposomes in nanostructured films, Mater. Sci. Eng., C, 2008.
40. Y. Y. Zuo, S. M. Tadayyon and E. Keating, et al., Atomic force microscopy studies of pulmonary surfactant films, Biophys. J., 2008, 2779–2791.
41. D. A. Doshi, A. M. Dattelbaum and E. B. Watkins, et al., Neutron reflectivity study of lipid membranes, Langmuir, 2005, 2865–2870.
42. J. D. Unsay, K. Cosentino and A. J. García-Sáez, Atomic force microscopy imaging and force spectroscopy, J. Visualized Exp., 2015, e52867.
43. T. J. Yacoub, A. S. Reddy and I. Szleifer, Structural effects and translocation of doxorubicin in lipid bilayers, Biophys. J, 2011, 378–385.
44. G. Minotti, C. Mancuso and A. Frustaci, et al., Cardiac lipid peroxidation in cancer patients treated with doxorubicin, J. Clin. Invest., 1996, 650–661.
45. H. Lee, Recent advances in simulation studies on the protein corona, Pharmaceutics, 2024, 1419.
46. J. Mariam, S. Sivakami and P. M. Dongre, Albumin corona on nanoparticles in drug delivery, Drug Delivery, 2016, 2668–2676.
47. T. Soranzo, D. K. Martin and J. L. Lenormand, et al., Coupling neutron reflectivity with cell-free protein synthesis, Sci. Rep., 2017, 3399.
48. D. Huster, H. A. Scheidt and K. Arnold, et al., Desmosterol may replace cholesterol in lipid membranes, Biophys. J., 2005, 1838–1844.
49. K. Roos, C. Wu and W. Damm, et al., OPLS3e: extending force field coverage, J. Chem. Theory Comput., 2019, 1863–1874.
50. L. Mohammed, H. Nourddine and E. F. Saad, et al., Chitosan-covered liposomes as drug transporter, RSC Adv., 2021, 1503–1516.

---