Pathway Complexity in Supramolecular Porphyrin Self-Assembly at an Immiscible Liquid–Liquid Interface

Nanostructures that are inaccessible through spontaneous thermodynamic processes may be formed by supramolecular self-assembly under kinetic control. In the past decade, the dynamics of pathway complexity in self-assembly have been elucidated through kinetic models based on aggregate growth by sequential monomer association and dissociation. Immiscible liquid–liquid interfaces are an attractive platform to develop well-ordered self-assembled nanostructures, unattainable in bulk solution, due to the templating interaction of the interface with adsorbed molecules. Here, we report time-resolved in situ UV–vis spectroscopic observations of the self-assembly of zinc(II) meso-tetrakis(4-carboxyphenyl)porphyrin (ZnTPPc) at an immiscible aqueous–organic interface. We show that the kinetically favored metastable J-type nanostructures form quickly, but then transform into stable thermodynamically favored H-type nanostructures. Numerical modeling revealed two parallel and competing cooperative pathways leading to the different porphyrin nanostructures. These insights demonstrate that pathway complexity is not unique to self-assembly processes in bulk solution and is equally valid for interfacial self-assembly. Subsequently, the interfacial electrostatic environment was tuned using a kosmotropic anion (citrate) in order to influence the pathway selection. At high concentrations, interfacial nanostructure formation was forced completely down the kinetically favored pathway, and only J-type nanostructures were obtained. Furthermore, we found by atomic force microscopy and scanning electron microscopy that the J- and H-type nanostructures obtained at low and high citric acid concentrations, respectively, are morphologically distinct, which illustrates the pathway-dependent material properties.


S5
S3. Model 2: coupling of competitive isodesmic and cooperative pathways. S7 S4. The Markov Chain Monte Carlo (MCMC) fitting procedure. S10 S5. Supplementary figures. S11 Figure S1. Experimental methodologies. S11 Figure S2. Baseline correction. S12 Figure S3. The influence of the aqueous phase pH on ZnTPPc interfacial selfassembly. S13 Figure S4. The influence of the porphyrin concentration on ZnTPPc interfacial selfassembly. S14 Figure S5. Scree plot and representation of the PCA results for the TIR-UV/vis dataset reported in Figure 1A (main text). S15 Figure S6. Multivariate Curve Resolution-Alternating Least Squares (MCR-ALS) analysis of the kinetics of interfacial ZnTPPc self-assembly. S16 Table S1. Quality control parameters extracted from MCR-ALS algorithms for each dataset. Table S2. Number of accepted and rejected steps, and number of iterations the covariance was updated for the MCMC run. S17 Figure S7. MCMC traces of the parameters from kinetic Model 1. S17 Figure S8. MCMC traces of the parameters from kinetic Model 2. S18 Figure S9. Histogram of the parameter distributions for kinetic Model 1. S19 Table S3. Parameter distributions obtained from the MCMC for kinetic Model 1. S20 Figure S10. Histogram of the parameter distributions for kinetic Model 2. S20 Table S4. Parameter distribution obtained from the MCMC for kinetic Model 2. S21 Table S5. Normalized sensitivity coefficients for kinetic Model 1. S22 Table S6. Normalized sensitivity coefficients for kinetic Model 2. S22 Figure S11. Sensitivity coefficients as a function of time for kinetic Model 1 for (A) the J-type nanostructure and (B) the H-type nanostructure. Figure S12. Sensitivity coefficients as a function of time for kinetic Model 2 for (A) the J-type nanostructure and (B) the H-type nanostructure. Figure S13. Scree plot and representation of the PCA results for the TIR-UV/vis dataset reported in Figure 6B (main text). In situ UV/vis spectroscopy in total internal reflection (TIR-UV/vis). ZnTPPc selfassembly was studied in situ at the interface between an aqueous phase containing 10 mM (analytical concentration) citric acid and a neat organic phase (TFT) by TIR-UV/vis using a custom-built optical setup ( Figure S1). Spectra were obtained every 0.5 s for up to 1000 s. The self-assembly process was performed using the "Porphyrin Last" protocol, 1  Ex situ microscopy characterisation. The porphyrin films were gently transferred to a silicon substrate for SEM and AFM analysis by bringing the solid support into contact with the interface. Prior to imaging, the samples were sequentially rinsed with water and TFT. AFM was performed using NT-MDT's Ntegra Spectra II. The topography was recorded using semicontact mode. The radius of curvature of the probe tip was less than 35 nm. The resonant frequency of the probe was 134.63 kHz. The probe stiffness was 5.83 N·m -1 . The gain of the lock-in-amplifier was set to 0.4. The scan size was set to 256 × 256 pixels. As the samples were expected to be rough in nature due to the stacking of flakes, the scan rate was 0.5 Hz.

Multivariate Curve Resolution-Alternating Least Squares (MCR-ALS) analysis.
MCR-ALS was used to analyze the spectral evolution of interfacial ZnTPPc at the aqueous|organic interface. This tool allows the spectral evolution of individual H-and J-type nanostructures to be separated. The calculations described herein were performed in R (the foundation for statistical computers, version 3.4.4). 3 All UV/vis spectra measured during the S4 self-assembly process were arranged in a data matrix Y(r × c), where rows r were the spectra recorded at different times during the reaction and columns c were kinetic profiles collected at different wavelengths. The raw spectra were further treated using the package baseline 4 for smoothing and correcting the drift of the signal ( Figure S2).
The decomposition of the matrix Y was achieved by using MCR-ALS, a wellestablished and widespread decomposition method meant to solve complex mixtures without any assumptions about the composition of the system. 5 Importantly, significant chemical information can be introduced in the optimization process under the form of constraints. 6,7 MCR-ALS decomposes the matrix Y according to the following equation where matrix S T (n × c) contains the spectral profiles of n pure resolved components, matrix C (r × n) describes the concentration profiles of these n species, and E represents the error matrix associated with the reconstruction. The first step of the decomposition required the determination of the number of significant components in the experimental data matrix by Principal Component Analysis (PCA). PCA is a reduction tool designed to identify the amount of variance present in a certain dataset and the determination of linearly independent components. Therefore, PCA allows the rank reduction of a large dataset to a few relevant components. This tool was implemented using the package FactoMineR 8 from R.
MCR-ALS was used according to the functions in package ALS. 9 For MCR, only the evolution of the porphyrin Soret band was analysed in the wavelength range from 400 to 484 nm. The initial matrix C was estimated by means of the detection of the purest concentration profiles 10 with a selected level of noise of 5%. The ALS routine was run employing the following soft constraints: non-negativity and unimodality for both pure spectra and concentration profiles. The "badness" of fit of the model obtained by MCR-ALS was evaluated by determining the Lack of Fit (LOF) parameter using the following equation: where and * are the experimental and calculated absorbance values, respectively. A second parameter was used to corroborate the quality of the optimization; the percentage of explained variance ( 2 ), defined as

S2. MODEL 1: COUPLING OF TWO COMPETITIVE COOPERATIVE (NUCLEATION-ELONGATION) PATHWAYS
The interfacial self-assembly of ZnTPPc can be explored using a one-dimensional nanostructure formation (or aggregation). This means we only consider the aggregation of a chain, and interactions with other chains in two-or three-dimensions are negligible. We modelled the size of the nanostructures by monomer association and dissociation, hence the proposed mechanism for the cooperative pathway (J-type nanostructure) is given by: The proposed mechanism for the other cooperative pathway (H-type nanostructure) is given by: To reduce the computational cost, the length of the polymers was set to 7. Further increasing the size of the chain did not affect the output. In the cooperative pathway, the last step is an autocatalytic process and therefore, to fulfil the mass balance, step 6 is irreversible following the model of Frieden for the polymerization of actin. 11 This mechanism produced a set of 13 ordinary differential equations (ODEs) that must be solved simultaneously. (S6) The initial conditions (numbers at t = 0) are as follows: The following constraints are applied: The mass balance of the reaction is: (S10) The corresponding system of coupled ODEs is solved using the functions provided by the deSolve package in R. 12 An interactive version of this model can be found in the following link https://entropia88.shinyapps.io/Shiny/

PATHWAYS
In this model, the interfacial self-assembly of ZnTPPc is explained by means of isodesmic and cooperative pathways. The isodesmic pathway is given by: The cooperative pathway is given by: To reduce the computational cost, the length of the polymers was set to 7. Further increases of the size of the chain did not affect the output. In the cooperative pathway, the last S9 step is an autocatalytic process. Therefore, to fulfil the mass balance, step 6 is irreversible, following the model of Frieden for polymerization of actin. 11 This mechanism produced a set of 13 ODEs that must be solved simultaneously. (S13) The initial conditions (numbers at t = 0) are as follows: The following constraints are applied: The mass balance of the reaction is: (S17) The corresponding system of coupled ODEs is solved using the functions provided by the deSolve package in R. 12 An interactive version of this model can be found at the following link https://entropia88.shinyapps.io/Shiny/

S4. THE MARKOV CHAIN MONTE CARLO (MCMC) FITTING PROCEDURE
The main issue affecting the resolution of bilinear data in MCR is the non-unicity of the solution, 13  A regular nonlinear fitting procedure is not the most accurate algorithm considering these feasible bands. Thus, it is important to provide an estimate of the parameter uncertainty, and to quantify the effect of that uncertainty on the observed variables. For this purpose, the S11 method chosen was Markov Chain Monte Carlo (MCMC). MCMC was implemented using the package FME. 12 FME uses MCMC with the Delayed Rejection and Adaptive Metropolis procedure. 15 MCMC is an efficient sampling method where the selection of the next parameter combination depends on the last parameter set and the resulting deviation between the model and the observation. Therefore, sampling is concentrated in the region with high likelihood.
This makes the method more efficient.
To avoid 'burn-in', the algorithm was started with the optimal parameter set as returned from the nonlinear fitting algorithm provided by the package FME. The parameter sets were taken using a Latin Hyperube Sampling method, with the Latinhyper function implemented in the FME package. MCMC was run with a number of 7000 steps, a delayed reaction of 2, the number of iterations after which the parameter covariance matrix is evaluated was set to 100, and a low variance weight was given to the prior distribution compared to the posterior distribution (wvar = 0.1). Additionally, the lower and upper bounds of the parameters were carefully selected to fulfil the conditions (S5) and (S12) in Models 1 and 2, respectively.  with the function baseline.rfbaseline for a Robust Baseline Estimation (a function that is included in the package baseline). 6 S13 ZnTPPc self-assembly at an immiscible aqueous|organic interface is very sensitive to the pH of the aqueous phase. 2,16 Above the pKa value of the carboxyl group on the porphyrin (pKa = 5.8), 9 the carboxyl groups are primarily deprotonated and the ZnTPPc monomers are negatively charged. Thus, the adsorbed monomers at the interface repel each other electrostatically and interfacial nanostructure self-assembly is inhibited. Consequently, only a band centred at 430 nm is detected in Figures S3D-F, corresponding to the Soret band of ZnTPPc monomers adsorbed at the interface. Below the pKa, the ZnTPPc molecules aggregate uncontrollably in the bulk aqueous phase (red lines, Figure S1A-B). Consequently, the quantity of free monomers adsorbed at the interface is very limited, and interfacial nanostructure formation does not take place. Only at a pH value equals to the pKa of the carboxyl groups can interfacial ZnTPPc self-assembly proceed ( Figure S3C). ZnTPPc at each pH. S14 Below 5 μM [ZnTPPc]aq., the interfacial self-assembly process did not proceed, and only the Soret band attributed to adsorbed monomers at the aqueous|organic interface were detected ( Figure S4A-B). For [ZnTPPc]aq. values between 5 and 10 μM, interfacial nanostructure formation proceeded, and all experimental datasets presented the same dynamic behaviour ( Figure S4C-E). This dynamic behaviour is discussed in detail in the main text and illustrated in Figure 1C. were identical (aqueous pH of 5.8, 10 mM citric acid aqueous electrolyte, and the organic phase was neat TFT). TIR-UV/vis spectra were taken every 0.5 s for 500 seconds. The red spectra are that of bulk aqueous ZnTPPc at pH 5.8. Figure S5A reports the scree plot for the TIR-UV/vis dataset in Figure 1A (main text; [ZnTPPc] of 4 nmol·cm -2 ). Eigenvalues associated to each component show an elbow in the second component, marking the separation between physiochemically meaningful components and noise-related ones. This result is confirmed by analyzing the scores and loadings from the PCA in Figure S5B-C. The scores and loadings appear noisier and unstructured from the third component. Based on these considerations, we ran the subsequent MCR-ALS analysis using two principal components (PCs, see Figure 2, main text). PC1 was associated to an H-type nanostructure (λmax. = 418 nm) and the PC2 to a J-type nanostructure (λmax. = 442 nm). The porphyrin monomers adsorbed at the interface were not a significant component that   Figure S6). Based on the pure spectra obtained, one spectrum was assigned to a J-type nanostructure (λmax. = 442 nm) and the other to a H-type nanostructure (λmax. = 418 nm). The peak that appeared at 452 nm corresponds to an artefact generated by small distortions of the spectra obtained in TIR Mode, and the small variations of the H-type nanostructure are attributed to the interference of the monomer spectrum (λmax 430 nm). Table S1 displays the quality parameters of the modelling, in all cases the % Lack of Fit (LOF) was close or below 5%, and the percentage of explained variance (r 2 ) was greater than 97%.

S19
The parameter distributions found by MCMC are plotted in Figures S9 and S10, and tabulated in Table S3 and S4, for Models 1 and 2, respectively. The dissociation constants of the nuclei for both models had a non-normal distribution, meaning that their values could change over a wide range, and this had no significant effect to the output. Hence, the uncertainty for these parameters was large. By contrast, the other parameters were normally distributed except for k1 and k3 in Model 2 which seem to have exponential distribution shapes.
The latter may be due to the lower bound restriction imposed on these parameters to fulfil condition (S12).   Figure S10. Histogram of the parameter distributions for kinetic Model 2.

S22
To further investigate the parameter uncertainty found by MCMC, a sensitivity analysis was done. Local sensitivity analysis was performed according to the function provided by the package FME. 12 A matrix, Sij, that contained the normalized and dimensionless sensitivities, was generated and every element of this matrix was defined as where yi is an output variable, and kj is the jth parameter. These coefficients were normalized by the nominal value of yi and kj. The higher the absolute sensitivity value, the more important the parameter. Thus, the magnitudes of the sensitivity can be used to rank the importance of the parameter to the output variable.
A summary of these ranks for Models 1and 2 are presented in Tables S5 and S6. Here, L1 and L2 are defined as follows:

S24
PCA was applied in order to determine the number of significant components in the dataset presented in Figure 6B (main text) where the aqueous citric acid concentration is set at 50 mM. From the scree plot in Figure S13A, the eigenvalues associated to the first component explain more than 97% of data variance. This result is further confirmed by analyzing the scores and loadings from the PCA ( Figure S13B-C). Although, from a simple inspection of Figure 6B the presence of more than one species is clearly discerned, the overlapping is so severe that a curve resolution will not successfully resolve the data. Based on these considerations, no MCR-ALS was performed on the datasets in Figure 6 involving large concentrations (50, 100 and 250 mM) of aqueous citric acid electrolyte.  Table S7. These results agree with the images obtained by SEM ( Figures 7A-B and 7E-F, main text).