Mechanism Study of Thermally Induced Anti-Tumor Drug Loading to Engineered Human Heavy-Chain Ferritin Nanocages Aided by Computational Analysis

Diverse drug loading approaches for human heavy-chain ferritin (HFn), a promising drug nanocarrier, have been established. However, anti-tumor drug loading ratio and protein carrier recovery yield are bottlenecks for future clinical application. Mechanisms behind drug loading have not been elaborated. In this work, a thermally induced drug loading approach was introduced to load anti-tumor drug doxorubicin hydrochloride (DOX) into HFn, and 2 functionalized HFns, HFn-PAS-RGDK, and HFn-PAS. Optimal conditions were obtained through orthogonal tests. All 3 HFn-based proteins achieved high protein recovery yield and drug loading ratio. Size exclusion chromatography (SEC) and transmission electron microscopy (TEM) results showed the majority of DOX loaded protein (protein/DOX) remained its nanocage conformation. Computational analysis, molecular docking followed by molecular dynamic (MD) simulation, revealed mechanisms of DOX loading and formation of by-product by investigating non-covalent interactions between DOX with HFn subunit and possible binding modes of DOX and HFn after drug loading. In in vitro tests, DOX in protein/DOX entered tumor cell nucleus and inhibited tumor cell growth.


Introduction
Mammalian ferritin is a 12 nm symmetrical protein cage consisting of 24 subunits. Each subunit contains a 4-helix bundle (helix A, B, C, and D) and a fifth short helix (helix E). Three N-terminals of subunits gather and form 8 hydrophilic channels in each ferritin shell to allow iron ion penetration [1]. Residues from 4 helices E make another 6 ferritin hydrophobic channels. All 14 channels on each ferritin shell are around 0.3-0.5 nm wide [2]. Ferritin's unique structure and high biocompatibility have made it a potential drug nanocarrier [3]. Especially, human heavy-chain ferritin (HFn) has shown an intrinsic active tumor targeting ability because it can recognize and bind to human transferrin receptor 1 (TfR1) [4].
Through decades of efforts, research have explored diverse drug loading approaches. Disassembly/reassembly and passive diffusion are 2 mainstream drug loading approaches. A disassembly/reassembly approach involves a dissociation/re-association of HFn assembly induced by pH or 8 M urea. This approach suits drugs either smaller or larger than ferritin channels. However, the disassembly/reassembly process in pH-induced pathway has been criticized. Kim et al. has proven that the process damaged ferritin structure and proteins (protein/DOX) stability test was performed to check drug leakage profile during storage. For the first time, computational analysis, molecular docking followed by MD simulation, was adopted to analyze chemical interactions contributing to drug loading and aggregation in thermally induced DOX loading process. Finally, in vitro evaluations, intracellular distribution, and cytotoxicity assays, compared 3 HFn-based proteins in vitro performances after thermally induced drug loading.

Materials
Three HFn-based proteins, HFn, HFn-PAS, and HFn-PAS-RGDK were designed as in a previous work [20]. Escherichia coli (E. coli) BL21 (DE3) (Tiangen Biotech, Beijing, China) was the expression host. MDA-MB-231 cell line was purchased from Cellbank Australia (Sydney, NSW, Australia). Cell culture related reagents were purchased from Thermo Scientific (Massachusetts, MA, USA). All other chemicals of analytical grade except for Doxorubicin hydrochloride (DOX) (Dalian Meilun Biotechnology, Dalian, China), were bought from Chem-Supply (Gillman, SA, Australia). All chromatography columns used in this work were bought from GE healthcare (Waukesha, WI, USA). Millipore purification system (Merck, Melbourne, VIC, Australia) was used throughout the experiments.

Preparation and Characterization of HFn and Functionalized HFns
E. coli strains expressing HFn or functionalized HFns were fermented in LB medium at 37 • C and target proteins were expressed by 0.5 mM isopropyl β-D-thiogalactoside (IPTG) 4 h induction. Harvested cell pellets were re-suspended, subjected to ultra-sonication for cell disruption. Lysis supernatants were collected and stored at −20 • C before purification. HFn was purified through the procedure established in a previous work [21]. Harvested E. coli lysis supernatants containing HFn-PAS and HFn-PAS-RGDK first underwent 50 • C, pH 5.0, 5 min heat-acidic precipitation to remove host cell proteins, buffer exchange using Hiprep X26/10 G25 desalting column (GE Healthcare, Waukesha, WI, USA), and then pH 7.0 mono Q ion-exchange chromatography (GE Healthcare, Waukesha, WI, USA) for polishing. The 12% reducing SDS-PAGE and TEM were adopted for purity and conformation integrity characterization, respectively. In TEM analysis, a FEI Tecnai G2 Spirit TEM (Eindhoven, NB, The Netherlands) was employed. Operating voltage was 100 kV. Three purified proteins were diluted to 0.1 mg mL −1 , spread on TEM support grids, air dried, and then negatively stained with 2% uranyl acetate before micrography capture.

Thermally-Induced Passive Loading of DOX into HFn, HFn-PAS-RGDK, and HFn-PAS
DOX was loaded to HFn-based nanocages through thermally induced passive diffusion. Temperature, buffer pH and incubation time are the main factors affecting drug loading. An orthogonal test was designed to optimize thermally induced drug loading condition for HFn and HFn-PAS-RGDK. Variables and levels tested are listed in Table 1. Initial protein concentration (1 mg mL −1 ) and DOX concentration (0.2 mg mL −1 ) were used in all conditions. Sample buffer was 20 mM phosphate buffer (PB) with 5 mM guanidinium chloride, pH 7.0 or 7.5. After thermal incubation of DOX and HFn-based nanocages, samples were at 1000 rpm 10 min at 4 • C to remove precipitates. Concentrations of the supernatants after centrifugation were measured using Bradford assay (Bio-Rad, Gladesville, NSW, Australia) for calculation of protein recoveries yields. Unloaded DOX was removed using Hitrap G25 desalting column (GE healthcare, Waukesha, WI, USA) and DOX loaded HFn-based protein (protein/DOX) were collected. All protein/DOX peaks then underwent SEC by Superose 6 increase 10/300 GL column (GE Healthcare, Waukesha, WI, USA) to detect if any soluble HFn-DOX aggregates existed.
SEC can separate DOX loaded in HFn-based nanocages (DOX loaded in nanocage) from soluble HFn-DOX aggregates. Peak areas (absorbance at 480 nm) can be used to determine the proportion of DOX loaded in nanocage, using Equation (1). Drug loading ratio, protein recovery yield, and the proportion of DOX loaded in nanocages under various conditions were compared to find the optimal condition. For HFn-PAS, DOX loading was conducted at the optimal loading condition of HFn-PAS-RGDK.

TEM Characterization of DOX Loaded HFn-Based Proteins and HFn-DOX Aggregate
Three protein/DOX samples under the optimal thermally induced drug loading conditions were analyzed using TEM analysis. A HFn-DOX aggregate sample collected from Superose 6 increase SEC also underwent TEM analysis. Sample treatment and device setting in TEM analysis were the same as in Section 2.2.

Computational Study of Interactions of HFn and DOX in Thermally-Induced Drug Loading
Molecular docking and Gromacs MD simulation analysis were used to identify the potential HFn and DOX interactions to explain the formation of HFn/DOX and soluble HFn-DOX aggregates. Molecular docking was performed to analyze the possible poses of HFn subunit and DOX interactions. MD simulation of the docking complexes aimed to find out the most stable HFn subunit-DOX complex structures.
Computational analysis was based on 2 prerequisites: (1) we assume that HFn subunit can be a representative of HFn assembly. This is because the assembly was theoretically 24 repetitions of the subunit. DOX is smaller than HFn channels, which makes it unlikely to simultaneously interact with more than one subunit of the same HFn assembly.
(2) The computational analysis focus was on the interactions between DOX and the residues Biosensors 2021, 11, 444 5 of 17 located on HFn assembly outer surface and inner surface. All interactions with interface residues of HFn assembly were ignored.
In molecular docking analysis, PyRx software was used. DOX 3D structure was from Pub Chem and HFn subunit structure file (PDB file) from RCSB PDB (ID: 2FHA). DOX was energy minimized before conducting docking. Top 9 docking HFn-DOX complexes were obtained and saved as PDB files. PDB files from docking results underwent Gromacs MD simulation using Gromacs 2018.
In MD simulation, CHARMM36 force field was used. The HFn-DOX complex structure was solvated in a dodecahedral box of size 460.73 nm 3 with water molecules and the box was charge neutralized by replacing eight water molecules with 8 Na+ ions. Energy minimization was conducted using the steepest descent integrator for 50,000 steps, until a tolerance of 10 kJ mol −1 . After this, temperature (NVT) and pressure equilibration (NPT) of the full system were performed at 323 K (approximate 50 • C). Finally, 10 ns 323 K simulation were conducted with 5,000,000 steps and 2 fs each step. Lincs constraint algorithm, Verlet cut-off scheme, Particle Mesh Ewald coulomb type were used in this MD simulation. Rootmean-square deviation (RMSD) and short-range non-bonded interaction energy of each complex in MD simulation were analyzed for the stability assessment. Three-dimensional structures of 9 complexes after MD simulation were saved and the interactions of HFn subunit and DOX within were visualized by Discovery Studio Visualizer. Interactions analyzed include hydrogen bond, salt bridge, and Pi (π) effects. Possible hydrophobic interaction was evaluated by analyzing the residue hydrophobicity in DOX binding area.

In Vitro Anti-Tumor Assessments of DOX Loaded HFn-Based Proteins
MDA-MB-231 is a human breast tumor cell line and has been proven to overexpress human TfR1, neuropilin 1, and integrin αvβ3/5 [22,23]. MDA-MB-231 cells were cultured in L-15 medium with 10% FBS and 1% PS. Intracellular distribution and MTT assay of all 3 protein/DOX were conducted.
MDA-MB-231 cells in exponential growth phase were utilized in intracellular distribution analysis and cytotoxicity assay. Procedures of these two assays were the same as in a previous work using another tumor cell line [20]. A well and cell viability were calculated using the following equations. IC 50 values of DOX and three protein/DOX were calculated in Origin 9.0 software. Unpaired T test was employed for statistical assessment. Figure 1 shows the SDS-PAGE and TEM images of 3 purified HFn-based proteins. In Figure 1A of 12% reducing SDS-PAGE, HFn subunit showed a single band with around 21 kDa. However, 2 functionalized HFns (HFn-PAS and HFn-PAS-RGDK) showed higher apparent molecular weights in SDS-PAGE gel than their theoretical 26 kDa and 26.5 kDa. The discrepancies in molecular weights probably result from PAS peptides, which has the tendency of binding to surrounding water molecules to increase the hydrodynamic radius. Other researchers discovered similar molecular weight increase in PAS modified proteins and PEG-conjugated proteins in SDS-PAGE analysis [18,24]. tendency of binding to surrounding water molecules to increase the hydrodynamic dius. Other researchers discovered similar molecular weight increase in PAS modi proteins and PEG-conjugated proteins in SDS-PAGE analysis [18,24]. TEM images in Figure 1B-D demonstrate that both functionalized HFns were ass bled hollow spheres, same as HFn. Cages of all proteins were around 12 nm in diam regardless of functionalization. This is because the inserted functional peptides at th terminus did not constitute the ferritin nanocage, while under TEM, the size of nanocage was visualized.

Optimization of HFn Thermally Induced Passive Loading to Increase Drug Loading
Thermally induced strategy takes advantage of the thermal energy mediated st tural perturbation of selective hydrophilic pore areas. Theil and co-workers used Circ Dichroism to analyze the α-helix content change of HFn following heat treatment at ferent temperatures, and found that a small amount of secondary structure began to t sition into random coil when temperature is greater than 45 °C, and it is very likely to place in pore areas and expand pores [25]. Heating also accelerates Brownian motio proteins and drug molecules so that greater efficiency could be achieved than in n heated passive diffusion. In this work, pH 7.0 and 7.5 were chosen to ensure that D carries positive charge (DOX pKa 8.3) and HFn inner surface has the opposite charge (H pI 4.8). Temperature conditions were selected based on thermal stability of HFn.
Standard curves for determination of drug loading ratio (N) are in Figure S1. Fig  2 summarizes the changes of drug loading ratios, proportions of DOX loaded in nanoc and HFn recovery yields with varying thermal induction time, temperature, and bu pH. Table S1 lists all drug loading ratios (Ns), proportions of DOX loaded into nanoc and protein recovery yields for HFn. As is shown in Figure 2A,C,E, with the increas thermal induction duration from 2 to 6 h, N increased at all tested temperatures. At 45 50 °C, and 60 °C, the highest N was 30.3, 41.6, and 56.7. Ns at pH 7.5 were slightly hig than those at pH 7.0 in most of the time regardless of temperature. TEM images in Figure 1B-D demonstrate that both functionalized HFns were assembled hollow spheres, same as HFn. Cages of all proteins were around 12 nm in diameter regardless of functionalization. This is because the inserted functional peptides at the C-terminus did not constitute the ferritin nanocage, while under TEM, the size of the nanocage was visualized.

Optimization of HFn Thermally Induced Passive Loading to Increase Drug Loading
Thermally induced strategy takes advantage of the thermal energy mediated structural perturbation of selective hydrophilic pore areas. Theil and co-workers used Circular Dichroism to analyze the α-helix content change of HFn following heat treatment at different temperatures, and found that a small amount of secondary structure began to transition into random coil when temperature is greater than 45 • C, and it is very likely to take place in pore areas and expand pores [25]. Heating also accelerates Brownian motion of proteins and drug molecules so that greater efficiency could be achieved than in non-heated passive diffusion. In this work, pH 7.0 and 7.5 were chosen to ensure that DOX carries positive charge (DOX pKa 8.3) and HFn inner surface has the opposite charge (HFn pI 4.8). Temperature conditions were selected based on thermal stability of HFn.
Standard curves for determination of drug loading ratio (N) are in Figure S1. Figure 2 summarizes the changes of drug loading ratios, proportions of DOX loaded in nanocage and HFn recovery yields with varying thermal induction time, temperature, and buffer pH. Table S1 lists all drug loading ratios (Ns), proportions of DOX loaded into nanocage and protein recovery yields for HFn. As is shown in Figure 2A,C,E, with the increase in thermal induction duration from 2 to 6 h, N increased at all tested temperatures. At 45 • C, 50 • C, and 60 • C, the highest N was 30.3, 41.6, and 56.7. Ns at pH 7.5 were slightly higher than those at pH 7.0 in most of the time regardless of temperature. In terms of proportions of DOX loading in nanocage, at 45 and 50 °C, they were at least 85% whilst at 60 °C they were below 85%. At all 3 temperatures, the proportion of DOX loaded into nanocage decreased with the duration of thermal induction. The proportions were largely pH-dependent, and the extent of pH-dependency was subject to temperature, hence they decreased by 0.4-3.5% at 45 °C and 50 °C, and by 10-15% at 60 °C as pH increased from 7.0 to 7.5.
For the HFn recovery yields, in Figure 2B,D,F, at 45 °C and 50 °C, they were above 90%, and at 60 °C, they were mostly below 85%. These results suggest that in the thermally induced drug loading process, DOX loaded in individual HFn nanocages, soluble HFn-DOX aggregates, and HFn-DOX precipitates were simultaneously produced as in previous research using disassembly/reassembly drug loading approaches. At 45 °C, proportion of drug loaded in nanocage and HFn recovery yield decreased slowly, and N increased slowly over time, whilst at 60 °C, proportion of drug loaded in nanocage, HFn recovery yield and N behaved in the opposite manner. These results confirm that 45 °C may not be effective to accelerate drug loading. In addition, at 60 °C, the local structures of HFn nanocages undergo excessive changes, resulting in massive formation of aggregates of HFn with DOX. Considering N, proportion of DOX loaded in nanocage and HFn recovery yield together, 50 °C, pH 7.5, 6 h is the best drug loading condition (N of 41.6, proportion of DOX loaded in nanocage of 87.2% and HFn recovery yield of 97.2%). In terms of proportions of DOX loading in nanocage, at 45 and 50 • C, they were at least 85% whilst at 60 • C they were below 85%. At all 3 temperatures, the proportion of DOX loaded into nanocage decreased with the duration of thermal induction. The proportions were largely pH-dependent, and the extent of pH-dependency was subject to temperature, hence they decreased by 0.4-3.5% at 45 • C and 50 • C, and by 10-15% at 60 • C as pH increased from 7.0 to 7.5.
For the HFn recovery yields, in Figure 2B,D,F, at 45 • C and 50 • C, they were above 90%, and at 60 • C, they were mostly below 85%. These results suggest that in the thermally induced drug loading process, DOX loaded in individual HFn nanocages, soluble HFn-DOX aggregates, and HFn-DOX precipitates were simultaneously produced as in previous research using disassembly/reassembly drug loading approaches. At 45 • C, proportion of drug loaded in nanocage and HFn recovery yield decreased slowly, and N increased slowly over time, whilst at 60 • C, proportion of drug loaded in nanocage, HFn recovery yield and N behaved in the opposite manner. These results confirm that 45 • C may not be effective to accelerate drug loading. In addition, at 60 • C, the local structures of HFn nanocages undergo excessive changes, resulting in massive formation of aggregates of HFn with DOX. Considering N, proportion of DOX loaded in nanocage and HFn recovery yield together, 50 • C, pH 7.5, 6 h is the best drug loading condition (N of 41.6, proportion of DOX loaded in nanocage of 87.2% and HFn recovery yield of 97.2%). Figure 3 and Table S2 show the DOX loading optimization results of HFn-PAS-RGDK. The relations between drug loading performance indicators (N, proportion of DOX loaded in nanocage and HFn-PAS-RGDK recovery yield) and experimental variables (induction time, pH and temperature) are similar to those in HFn.

Optimization of HFn-PAS-RGDK Thermally Induced Passive Loading
Biosensors 2021, 11, x FOR PEER REVIEW 8 of 17 Figure 3 and Table S2 show the DOX loading optimization results of HFn-PAS-RGDK. The relations between drug loading performance indicators (N, proportion of DOX loaded in nanocage and HFn-PAS-RGDK recovery yield) and experimental variables (induction time, pH and temperature) are similar to those in HFn. In Figure 3A,C,E,F positively responded to thermal induction duration and temperature. pH 7.5 showed greater Ns than pH 7.0 in most of time. Proportion of DOX loaded in nanocage was negatively related to temperature and incubation time. At 45 °C and 50 °C , proportions of DOX loaded in nanocage were greater than 75%. At 60 °C, they were lower than 70%. As in Figure 3B,D,E, HFn-PAS-RGDK recovery yields were greater than 75% except at 60 °C 4 h and 6 h. The best DOX loading condition was obtained at 50 °C, pH 7.5 and 6 h, with an N of 45.2, proportion of DOX loaded in nanocage of 78.5% and HFn-PAS-RGDK recovery yield of 76.0%. HFn-PAS DOX loading ratio was 38.4, proportion of DOX loaded in nanocage was 73.4% and protein recovery was 75.1% under the same condition.

Conformation of DOX Loaded HFn and Functionalized HFns
SEC analysis was performed to prove the success of DOX loading and separate DOX loaded nanocages from HFn-DOX soluble aggregates according to hydrodynamic volume differences. HFn-based proteins have absorbance at 280 nm but not at 480 nm. DOX has In Figure 3A,C,E, N positively responded to thermal induction duration and temperature. pH 7.5 showed greater Ns than pH 7.0 in most of time. Proportion of DOX loaded in nanocage was negatively related to temperature and incubation time. At 45 • C and 50 • C, proportions of DOX loaded in nanocage were greater than 75%. At 60 • C, they were lower than 70%. As in Figure 3B,D,E, HFn-PAS-RGDK recovery yields were greater than 75% except at 60 • C 4 h and 6 h. The best DOX loading condition was obtained at 50 • C, pH 7.5 and 6 h, with an N of 45.2, proportion of DOX loaded in nanocage of 78.5% and HFn-PAS-RGDK recovery yield of 76.0%. HFn-PAS DOX loading ratio was 38.4, proportion of DOX loaded in nanocage was 73.4% and protein recovery was 75.1% under the same condition.

Conformation of DOX Loaded HFn and Functionalized HFns
SEC analysis was performed to prove the success of DOX loading and separate DOX loaded nanocages from HFn-DOX soluble aggregates according to hydrodynamic volume differences. HFn-based proteins have absorbance at 280 nm but not at 480 nm. DOX has absorbance at both wavelengths. Protein/DOX in theory has absorbance at both wavelengths and peak retention time should be similar to HFn-based proteins. Figure 4A-C show the chromatograms of HFn/DOX, HFn-PAS/DOX, and HFn-PAS-RGDK/DOX prepared under the condition of 50 • C, 6 h, pH 7.5. Two peaks, P1 and P2, were observed in all 3 samples. The larger P2 had retention volumes of 13-15 mL in Superose 6 increase column, and absorbance at both 280 and 480 nm. This means it was the DOX loaded HFn-based nanocage. The smaller P1 eluted before P2 was the protein-DOX soluble aggregates, the DOX amount of which accounted for below 27% of the total DOX in the SEC loading samples. The heating process did not affect most of the ferritin nanocage, as are shown in Figure 4D-F. Most of the protein/DOX were hollow spheres. Nanocage sizes were still around 12 nm, the same as before thermally induced passive drug loading process.
Biosensors 2021, 11, x FOR PEER REVIEW 9 of 17 absorbance at both wavelengths. Protein/DOX in theory has absorbance at both wavelengths and peak retention time should be similar to HFn-based proteins. Figure 4A-C show the chromatograms of HFn/DOX, HFn-PAS/DOX, and HFn-PAS-RGDK/DOX prepared under the condition of 50 °C, 6 h, pH 7.5. Two peaks, P1 and P2, were observed in all 3 samples. The larger P2 had retention volumes of 13-15 mL in Superose 6 increase column, and absorbance at both 280 and 480 nm. This means it was the DOX loaded HFnbased nanocage. The smaller P1 eluted before P2 was the protein-DOX soluble aggregates, the DOX amount of which accounted for below 27% of the total DOX in the SEC loading samples. The heating process did not affect most of the ferritin nanocage, as are shown in Figure 4D-F. Most of the protein/DOX were hollow spheres. Nanocage sizes were still around 12 nm, the same as before thermally induced passive drug loading process.

DOX Loaded HFn and Functionalized HFns Stability
Protein/DOX stability test was designed to reflect the stability of protein/DOX in storage (4 °C) and blood circulation (37 °C). In storage, drug leakage profiles for all protein/DOX were consistent, where around 20% of loaded drug leaked over 2 weeks ( Figure  5). At 37 °C, protein/DOX were less stable in contrast with 4 °C, with around 30% of drug loss detected in HFn/DOX and 35% of drug leaking in other 2 groups for 2 weeks. In all 3 protein/DOX, drug leaked fast during the initial 12 h, then slowed down. Perhaps some of the loaded drugs were just loosely attached or physically trapped inside protein nanocages. Hence, these drugs were more prone to dissociation from protein, while drugs strongly interacted with HFn remained within ferritin. Functionalized HFns showed lower protein/DOX stabilities compared with HFn, which probably result from the insertion of foreign peptides.

DOX Loaded HFn and Functionalized HFns Stability
Protein/DOX stability test was designed to reflect the stability of protein/DOX in storage (4 • C) and blood circulation (37 • C). In storage, drug leakage profiles for all protein/DOX were consistent, where around 20% of loaded drug leaked over 2 weeks (Figure 5). At 37 • C, protein/DOX were less stable in contrast with 4 • C, with around 30% of drug loss detected in HFn/DOX and 35% of drug leaking in other 2 groups for 2 weeks. In all 3 protein/DOX, drug leaked fast during the initial 12 h, then slowed down. Perhaps some of the loaded drugs were just loosely attached or physically trapped inside protein nanocages. Hence, these drugs were more prone to dissociation from protein, while drugs strongly interacted with HFn remained within ferritin. Functionalized HFns showed lower protein/DOX stabilities compared with HFn, which probably result from the insertion of foreign peptides. loss detected in HFn/DOX and 35% of drug leaking in other 2 groups for 2 weeks. In all 3 protein/DOX, drug leaked fast during the initial 12 h, then slowed down. Perhaps some of the loaded drugs were just loosely attached or physically trapped inside protein nanocages. Hence, these drugs were more prone to dissociation from protein, while drugs strongly interacted with HFn remained within ferritin. Functionalized HFns showed lower protein/DOX stabilities compared with HFn, which probably result from the insertion of foreign peptides.

Interactions between HFn and DOX in Thermally Induced Drug Loading by Computational Analysis
From molecular docking results, 9 different HFn subunit-DOX complexes were obtained. Complexes underwent 10 ns 50 • C MD simulation for stability assessment. Threedimensional structures of 9 complexes after simulation are shown in Figure 6. Among them, only in Complex 1, the location DOX binds to was the inner surface in HFn assembly. In the other 8 complexes, DOX bound to areas corresponding to the outer surface in HFn assembly. This implies that Complex 1 is very likely to be the structure of DOX loaded in HFn nanocage, while the interaction ways in the other complexes could form drug loading, soluble HFn-DOX aggregates, and HFn-DOX precipitates in thermally induced drug loading process.

Interactions between HFn and DOX in Thermally Induced Drug Loading by Computational Analysis
From molecular docking results, 9 different HFn subunit-DOX complexes were obtained. Complexes underwent 10 ns 50 °C MD simulation for stability assessment. Threedimensional structures of 9 complexes after simulation are shown in Figure 6. Among them, only in Complex 1, the location DOX binds to was the inner surface in HFn assembly. In the other 8 complexes, DOX bound to areas corresponding to the outer surface in HFn assembly. This implies that Complex 1 is very likely to be the structure of DOX loaded in HFn nanocage, while the interaction ways in the other complexes could form drug loading, soluble HFn-DOX aggregates, and HFn-DOX precipitates in thermally induced drug loading process. To evaluate the stabilities of these structures, the RMSD and the short-range nonbonded interaction energy of HFn subunit and DOX molecule in 9 complexes during simulation were monitored and are presented in Figures 7 and 8. The smaller the RMSD and the lower the energy is, the more stable the complex structure is and the more reliable the complex structure is. The stability orders of structures demonstrated in RMSD and energy are consistent. To evaluate the stabilities of these structures, the RMSD and the short-range nonbonded interaction energy of HFn subunit and DOX molecule in 9 complexes during simulation were monitored and are presented in Figures 7 and 8. The smaller the RMSD and the lower the energy is, the more stable the complex structure is and the more reliable the complex structure is. The stability orders of structures demonstrated in RMSD and energy are consistent.     Complex 1 was the most stable structure in 50 • C MD simulation, with RMSD lower than 1 nm and energy below −350 kJ mol −1 . This result is in accordance with the experiment result conducted under 50 • C, that more DOX was being loaded to HFn nanocage than forming soluble aggregates and precipitates. Complex 4 and 5 were the second most stable structures, of which the RMSD were below 1 nm and the energies were below −200 kJ mol −1 most of the time. Complex 3 was the third most stable structure. Its RMSD was lower than 1.2 nm. RMSD of Complex 7,6,8,9, and 2 were greater than 1.5 nm and energies of them were above −50 kJ mol −1 , indicating relatively unstable structures and possibly weak interactions. Based on the stability results, Complex 1, 4, 5, and 3 were the focus in interaction analysis. Table 2 lists the hydrogen bond, salt bridge and Pi effect interactions between HFn subunit and DOX in Complex 1, 4, 5, and 3 at 10 ns of the simulation. Figure S4 lists the 2D diagrams of 9 complexes of HFn subunit with DOX at 10 ns of the MD simulation. Hydrogen bonds and salt bridges are strong non-covalent bonds, in contrast with van der Waals interaction, such as Pi effects. The more of them in the complex, the more stable the complex structure is. Complex 1 had the most hydrogen bonds and salt bridges.  1 Numbers in the brackets after the residue are the number of the interactions involving the residue.
Regarding the possible hydrophobic interactions between HFn subunit and DOX in Complex 1, 4, 5, and 3, a 5 residue average hydrophobicity was used to reflect the hydrophobicity level of residues in DOX binding area. This is because in a protein, the hydrophobicity of residues can be affected by the nearby residues. Local area hydrophobicity reflects the possibly of hydrophobic interaction better than considering individual residue hydrophobicity. The calculation of 5 residue average hydrophobicity has considered the impact of nearby residues and its value demonstrates how hydrophobic the local area of the residue is. Hydrophobicity values in Table 3 were calculated using Discovery Studio Visualizer. The greater the value is, the more likely it would interact with the hydrophobic DOX molecule. In Complex 1, 4, 5, and 3, there were at least 3 residues at the binding pocket available for hydrophobic interactions with DOX. Table 3. Five residue average hydrophobicity of residues in DOX binding area in complex 1, 4, 5, and 3. According to the computational analysis, in DOX loaded HFn nanocage, DOX was mostly bound to HFn subunit as in Complex 1. Relatively weak binding ways found in Complex 4, 5, and 3 and physically trapped DOX also existed. Therefore, the loading ratio could reach above 24. However, DOX remained in HFn in these weaker ways are more prone to dissociation during storage. Physically trapped DOX probably accounts for the burst release of DOX in the initial 12 h, and the weakly bounded DOX on HFn surface in Complex 4, 5, and 3 would gradually be released, as is shown in Figure 5.

Complex 5 Residue Average Hydrophobicity Values of Hydrophobic Residues in DOX
Because no aggregates nor precipitates were found in 50 • C 6 h heated HFn, it is reasonable to infer that the interaction of DOX and HFn assembly has led to HFn and DOX aggregation. Small aggregates are still soluble while huge ones turn into HFn-DOX precipitates. TEM image in Figure 9A demonstrates that the HFn in soluble HFn-DOX aggregates were still intact spheres but clumped into a large particle. Interaction ways in Complex 4, 5, and 3 and others, except Complex 1, are theoretically possible to cause aggregation in a way that DOX works as a cross linker ( Figure 9B). In each HFn assembly, there are 24 subunits for DOX to bind to, and the HFn-DOX aggregates contain multiple DOX molecules.
According to the computational analysis, in DOX loaded HFn nanocage, DOX w mostly bound to HFn subunit as in Complex 1. Relatively weak binding ways found Complex 4, 5, and 3 and physically trapped DOX also existed. Therefore, the loading ra could reach above 24. However, DOX remained in HFn in these weaker ways are m prone to dissociation during storage. Physically trapped DOX probably accounts for burst release of DOX in the initial 12 h, and the weakly bounded DOX on HFn surface Complex 4, 5, and 3 would gradually be released, as is shown in Figure 5.
Because no aggregates nor precipitates were found in 50 °C 6 h heated HFn, it is r sonable to infer that the interaction of DOX and HFn assembly has led to HFn and D aggregation. Small aggregates are still soluble while huge ones turn into HFn-DOX p cipitates. TEM image in Figure 9A demonstrates that the HFn in soluble HFn-DOX agg gates were still intact spheres but clumped into a large particle. Interaction ways in Co plex 4, 5, and 3 and others, except Complex 1, are theoretically possible to cause aggre tion in a way that DOX works as a cross linker ( Figure 9B). In each HFn assembly, th are 24 subunits for DOX to bind to, and the HFn-DOX aggregates contain multiple D molecules.

Intracellular Distribution and Cytotoxicity of DOX Loaded HFn-Based Proteins
DOX has been proven to be able to diffuse into cell nucleus and disrupt cell divis [26]. However, in theory, the DOX loaded on HFn and functionalized HFns need to released from protein prior to exerting its function. Intracellular distribution test aimed check whether the release of DOX form protein/DOX occurred. Cell nucleus locatio were visualized as blue dots under cell imager after Hoechst 33,258 staining ( Figure 10 Due to the intrinsic fluorescence of DOX molecules, under the excitation of 480 nm lig DOX molecule accumulation could be observed as green dots. In the merged images, color of dots in all four groups changed to light cyan, indicating that DOX molecu loaded on proteins through thermally induced loading approach were released and ac mulated inside cell nucleus.

Intracellular Distribution and Cytotoxicity of DOX Loaded HFn-Based Proteins
DOX has been proven to be able to diffuse into cell nucleus and disrupt cell division [26]. However, in theory, the DOX loaded on HFn and functionalized HFns need to be released from protein prior to exerting its function. Intracellular distribution test aimed to check whether the release of DOX form protein/DOX occurred. Cell nucleus locations were visualized as blue dots under cell imager after Hoechst 33,258 staining ( Figure 10A). Due to the intrinsic fluorescence of DOX molecules, under the excitation of 480 nm light, DOX molecule accumulation could be observed as green dots. In the merged images, the color of dots in all four groups changed to light cyan, indicating that DOX molecules loaded on proteins through thermally induced loading approach were released and accumulated inside cell nucleus.
MTT assay was designed to compare the inhibition effects of DOX loaded on HFn, two functionalized HFns, and free DOX. Figure 10B shows the cell viabilities at different concentrations of equivalent DOX. Table 4 lists the IC 50 values of all four groups. Free DOX possessed the lowest IC 50 . However, it does not indicate free DOX has the greatest anti-proliferation effect. This is because in in vitro assays, the direct incubation of free DOX with cells has maximized the internalization efficiency of free DOX. On the contrary, the uptake efficiency of DOX in in vivo tests and in real practice would be greatly hampered by the blood circulation and metabolism system. Except for free DOX, DOX loaded on HFn-PAS-RGDK had the lowest IC 50 value, followed byHFn-PAS/DOX. HFn/DOX had the greatest IC 50 . MTT assay was designed to compare the inhibition effects of DOX loaded on HFn, two functionalized HFns, and free DOX. Figure 10B shows the cell viabilities at different concentrations of equivalent DOX. Table 4 lists the IC50 values of all four groups. Free DOX possessed the lowest IC50. However, it does not indicate free DOX has the greatest anti-proliferation effect. This is because in in vitro assays, the direct incubation of free DOX with cells has maximized the internalization efficiency of free DOX. On the contrary, the uptake efficiency of DOX in in vivo tests and in real practice would be greatly hampered by the blood circulation and metabolism system. Except for free DOX, DOX loaded on HFn-PAS-RGDK had the lowest IC50 value, followed byHFn-PAS/DOX. HFn/DOX had the greatest IC50.
No statistical significance was found between anti-proliferation abilities of HFn-PAS-RGDK/DOX group and free DOX group (p > 0.05). Anti-proliferation effect of HFn-PAS-RGDK/DOX was significantly higher than the other two HFn-based protein/DOX groups (p < 0.05). This is because of the tumor targeting ability of the inserted RGDK in HFn-PAS-RGDK.

Discussion
In this study, the thermally induced passive diffusion approach succeeded in loading DOX into HFn and 2 functionalized HFns. 50 °C, pH 7.5 and 6 h was found to be the optimal condition for HFn and functionalized HFns. Temperature and incubation time showed a great impact on DOX loading performance. Although HFn and DOX have outstanding thermal stabilities, in the thermally induced drug loading process, both drug loading and irreversible HFn-DOX aggregation occurred under all selected conditions. With the same incubation time, as incubation temperature increased, N value increased whilst the proportion of DOX loaded in nanocage declined. At the same incubation temperature, N value increased and the proportion of DOX loaded in nanocage decreased over incubation time, especially at 60 °C.  No statistical significance was found between anti-proliferation abilities of HFn-PAS-RGDK/DOX group and free DOX group (p > 0.05). Anti-proliferation effect of HFn-PAS-RGDK/DOX was significantly higher than the other two HFn-based protein/DOX groups (p < 0.05). This is because of the tumor targeting ability of the inserted RGDK in HFn-PAS-RGDK.

Discussion
In this study, the thermally induced passive diffusion approach succeeded in loading DOX into HFn and 2 functionalized HFns. 50 • C, pH 7.5 and 6 h was found to be the optimal condition for HFn and functionalized HFns. Temperature and incubation time showed a great impact on DOX loading performance. Although HFn and DOX have outstanding thermal stabilities, in the thermally induced drug loading process, both drug loading and irreversible HFn-DOX aggregation occurred under all selected conditions. With the same incubation time, as incubation temperature increased, N value increased whilst the proportion of DOX loaded in nanocage declined. At the same incubation temperature, N value increased and the proportion of DOX loaded in nanocage decreased over incubation time, especially at 60 • C. Table 5 compares HFn drug loading performance of this work with some previously published studies. In this study, N of HFn (41.6) is greater than previously studies, which adopted 8 M urea or optimized stepwise pH-induced disassembly-reassembly approaches. Recovery yield of HFn in this study, 97.2%, is similar to high hydrostatic pressure passive diffusion approach (99%) and significantly greater than the pH-induced (25%, 55%) or 8 M urea-based approach (64.8%). Disassembly/reassembly approach has been questioned to be not fully reversible because 2 holes were detected by synchrotron small-angle X-ray scattering, and the authors argued that this structural damage may result in protein loss and aggregation in the drug loading process [5]. To the contrary, at 50 • C, HFn nanocage remains intact throughout the thermally induced drug loading process, which involves less structural changes. '/' means no data.
Compared with HFn, under most experimental conditions, especially at 50 • C and 60 • C, the functionalized HFn, HFn-PAS-RGDK, had relatively low protein recovery yields and low proportions of DOX loaded in nanocage. HFn-PAS also demonstrated reduced proportions of DOX loaded in nanocage and protein recovery yields. Two functionalized HFns were more prone to aggregation in the heating process, suggesting slightly decreased thermal stabilities. This could be ascribed to the 'flip to flop' phenomenon in functionalized HFns, where E-helix with inserted functional peptide are extruded outside HFn nanocage, as was discovered in our previous work [20]. Hydrophobic interactions of 4 helices E around each hydrophobic channel in natural 'flip' HFn have been proven to contribute to HFn stability [29,30]. The turnover of E-helix has hampered helices E interactions and thus declined thermal stability.
Combining the results from molecular docking, MD simulation, and experiments, hydrogen bond and salt bridges between DOX and HFn residues in Complex 1 probably account for most of the loading of DOX. Physical entrapment of DOX in HFn assembly and interactions in other complexes may also contribute but they suffer from a rapid DOX leakage during storage, as shown in Figure 5. In the process of thermally induced DOX loading, DOX may undergo unexpected interactions with multiple HFn assemblies through hydrogen bonds and salt bridges to form HFn-DOX aggregates ( Figure 9B).
In vitro tests demonstrate that DOX loaded through thermally induced passive diffusion could exert anti-cancer function as free DOX.

Conclusions
A thermally induced drug loading approach has improved DOX loading ratios and protein recovery yields of HFn and functionalized HFns, HFn-PAS and HFn-PAS-RGDK. This mild and efficient strategy can become an alternative to produce HFn-based nanocages with various drugs. According to molecular docking and MD simulation analysis, hydrogen bond, salt bridges and other non-covalent interactions between HFn and DOX molecules contribute to DOX loading and by-product formation. The combination of molecular docking and MD simulation analyses can be a useful tool to shed light on ferritin drug loading mechanism. In vitro tests show that thermally-induced DOX loaded HFn-based proteins can exert tumor inhibition of DOX. RGDK has promoted DOX internalization to tumor cells and enhanced HFn anti-tumor efficacy.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/bios11110444/s1, Figure S1: Standard linear curves of correlations between drug or HFn-based protein nanocages concentrations and optical densities. Table S1: Loading ratios (Ns), proportions of DOX loaded in nanocage and protein recovery percentages in HFn thermally induced drug loading optimization. Table S2: Loading ratios (Ns), proportions of DOX loaded in nanocage and protein recovery percentages in HFn-PAS-RGDK thermally induced drug loading optimization.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.