Turing patterns by supramolecular self-assembly of a single salphen building block.

In the 1950s, Alan Turing showed that concerted reactions and diffusion of activating and inhibiting chemical species can autonomously generate patterns without previous positional information, thus providing a chemical basis for morphogenesis in Nature. However, access to these patterns from only one molecular component that contained all the necessary information to execute agonistic and antagonistic signaling is so far an elusive goal, since two or more participants with different diffusivities are a must. Here, we report on a single-molecule system that generates Turing patterns arrested in the solid state, where supramolecular interactions are used instead of chemical reactions, whereas diffusional differences arise from heterogeneously populated self-assembled products. We employ a family of hydroxylated organic salphen building blocks based on a bis-Schiff-base scaffold with portions responsible for either activation or inhibition of assemblies at different hierarchies through purely supramolecular reactions, only depending upon the solvent dielectric constant and evaporation as fuel.

Introduction

In Turing’s original paper (Turing, 1952), the so-called 6th morphology consists of stationary periodic waves with finite wavelength (λ), where the crests correspond to the visible part in the pattern (spots, stripes, or labyrinths) and the troughs in the same are their respective void spaces. In Nature, Turing patterns (e.g. Figure 1) can be found in the colony growth and spatial positioning mechanisms of bacteria (Karig et al., 2018; Murray and Sourjik, 2017), pigment distribution and structural features in fish (Almuedo-Castillo et al., 2018; Cooper et al., 2018; Konow et al., 2021; Mahalwar et al., 2014; Onimaru et al., 2016), seashells (Cooper, 2012; Meinhardt, 2009), corals (Cohen et al., 2004; Pratchett et al., 2015), octopus (Ishida, 2021), plants (Vadde and Roeder, 2020; van den Berg et al., 2021), anatomical traits and hair follicle spacing in mammals (Cetera et al., 2018; Economou et al., 2012, 2020; Raspopovic et al., 2014; Sheth et al., 2012; Sick et al., 2006), feather branching and coloration in birds (Harris et al., 2005; Haupaix et al., 2018), the left-right asymmetry and teeth development in vertebrates (Marcon and Sharpe, 2012; Salazar-Ciudad and Jernvall, 2010), etc.

Figure 1

Turing patterns are ubiquitous in Nature

(A) The coral Diploria labyrinthiformis which is usually present in the Atlantic Ocean exhibits superficial 2D waves that meet Turing’s criteria; the specimen in the figure was found by the authors as a brick in a wall of the first house that Hernán Cortés constructed during the XVI century in La Antigua, Veracruz, central east Mexico; scale bar is 10 mm, inset is the full specimen in its current location at the coordinates 19°19′16.2"N, 96°19′15.1"W.

(B) Fast Fourier transform (FFT) of the photograph shows a 360° annular spot rendered by the isotropic wavelength in the pattern (λ), where the tolerance in λ is indicated by the rim of the ring’s width limits, i.e. for the case, 5.6–10.8 mm, while the central value in λ is represented by the circumference of maximum intensity located between these limits (λ = 8.2 ± 2.6 mm).

Turing pattern formation mechanism is based on, at least, two coexisting morphogens, an activator (A) and an inhibitor (B), where diffusion (DA, DB …) can lead to patterning if antagonistic reactions at the long-range are coupled to agonistic reactions at the short- to the mid-range (Gierer and Meinhardt, 1972; Kondo and Miura, 2010; Turing, 1952). Commonly, this means that a faster degradation (or a slower synthesis) of A is remotely driven by the larger mobility of B, whereas the morphogenetic output (the patterning readout) is locally triggered by the less mobile A counterpart; this usually implies that an easier patterning is possible for the larger the difference between DA and DB. Such an interplay between reactions and diffusion of A and B occurs through stochastic instabilities arisen from initially stable and homogeneous concentrations. In this regard, the regulatory role of the components is associated with non-linearity in the reactions and cooperativity, in a way that for the larger the cooperativity in a reaction-diffusion (RD) system, the lesser the restrictions for Turing patterns to emerge (Diambra et al., 2015; Epstein and Xu, 2016; Tompkins et al., 2014).

These principles have been observed not only in biological morphogenetic processes but are also widely present in many more dynamic systems than thought earlier (Kuznetsov and Polezhaev, 2020; Leyshon et al., 2021; Scholes et al., 2019), including demographic, sociolinguistic, psychologic, economic, ecologic, or epidemiologic phenomena (Batabyal, 2021; Chakraborty et al., 2021; Chen, 2019; Iskarous, 2019; Lacalli, 2020; Mimar et al., 2021; Pal and Poria, 2021; Putra et al., 2019; Vandermeer and Perfecto, 2020; Zincenko et al., 2021). Nevertheless, in chemistry, these RD principles have been typically used to explain how reactions concerted at different scales can result into oscillatory and out-of-equilibrium systems, such as the Belousov-Zhabotinsky (BZ) reaction (Zaikin and Zhabotinsky, 1970) (e.g. oxidation of malonic acid by bromate using metal-catalysis) and some other examples in biochemistry and metabolism, organic chemistry, and redox and electrochemical reactions (Agladze and Steinbock, 2000; Astrov and Logvin, 1997; Christoph et al., 1999; Ertl, 2008; Fleury, 1997; Goesten et al., 2016; He et al., 2012; Heuser et al., 2015; Nakouzi and Sultan, 2011; Sugai et al., 2017; Yashin and Balasz, 2006).

So far, the RD principles have been more widely employed to engineer synthetic systems based on (covalent) chemical reaction networks. These systems facilitated the access to autonomous evolutionary synthesis of functional cages, macrocycles and molecular motors, self-replication of abiotic molecules, and bottom-up construction of adaptive biosensors and self-assembled architectures (Cera and Schalley, 2018; Grzelczak et al., 2019; Grzybowski et al., 2017; Kosikova and Philp, 2017; Orrillo et al., 2019; Scalise and Schulman, 2013; van Roekel et al., 2015; van Rossum et al., 2017). However, it was not until 2021, that such principles were experimentally demonstrated to also allow for soft-/inorganic-matter redistribution patterning that otherwise were believed to be microphase separation or nucleation-growth processes, with the construction of out-of-the-equilibrium Turing patterns arrested in the solid state at the macro, the micro, and the nanoscale in polymeric membranes (Qiao et al., 2021; Wu et al., 2021), semiconducting wafers (Leitgeb et al., 2021), atomic layer deposits (Fuseya et al., 2021), and metal chalcogenide nanostructures (Zhang et al., 2021). These works proved the existence of an interplay between RD in Turing instabilities and classical phase-separation processes formerly modeled through the Cahn-Hilliard approach, which on the contrary describes processes occurring near the thermal equilibrium. In this manner, and as earlier predicted analytically (Brauns et al., 2020, 2021), the reports also afforded a wider perspective on RD patterning, showing that morphogenesis also proceeds through mass redistribution instabilities driven by multiple local equilibria that play as mass-flux scaffolds to which the pattern is intrinsically linked. Unfortunately, the use of Turing’s model to engineer patterning systems based on much less complex RD platforms beyond the reported examples or by using “supramolecular reactions” instead of “covalent reactions” as a paradigm shift (Vantomme and Meijer, 2019; Whitesides et al., 1995) that, in turn, could help in reconceiving the classical physical point of view for matter aggregation into a more inclusive frame that accounted for non-covalent patterning from small molecules has remained a major challenge over the last decades (Bánsági et al., 2011; Grzybowski et al., 2017; Kondo and Miura, 2010; Nakouzi and Steinbock, 2016; Santos-Moreno and Schaerli, 2019). In this sense, autonomous pattern evolution from a single molecular structure which itself contained all the information needed to cooperatively achieve auto-regulation via agonistic and antagonistic signaling and, consequently, tuneable self-organization at increased complexity levels would represent a breakthrough in matter control and adaptive chemistry (Lehn, 2013), since no reports exist to date in this regard, presumably by the challenge in programming antagonistic “covalent” reactivities in a same synthon. The latter would be a single-molecule level example of how equal cells can self-organize and spatially constrain into differentially controlled structures in living organisms without the need for external signals but relying only on collective information exchanges, as pointed out by the latest discoveries in biology (Marcon and Sharpe, 2012; Raspopovic et al., 2014; Sheth et al., 2012; Toda et al., 2018) yet artificially shown so far only for self-sorted DNA origami nanofibrils and nanotubes lacking spatial distribution control (Groeer et al., 2021). Therefore, instead of using an orthodox approach based on agonist/antagonist “classical” (covalent) chemical reactions to transfer information between molecules of different predetermined diffusivities, herein, we used supramolecular interactions as the signaling information carriers and differently aggregated populations to engineer an artificial RD system with autonomous organization capabilities through self-assembly of a single building block (BB). In this way, we simplified our task to the identification of the minimum molecular characteristics and medium conditions that could facilitate the access to Turing patterns arrested in the solid state as a bottom-up construction tool.

In a recent work (Zelada-Guillén et al., 2018), we studied the tuneable self-assembly of modular BBs based on bis-salphen compounds, i.e. symmetrical tetrakis-Schiff bases containing two terminal di-aromatic ketimine moieties linked through adjacent o-phenylenediamine derivative portions, both being coupled via aldimine counterparts to a central biphenyl hinge. In that report, we demonstrated that those structures could drive concerted self-assembly scenarios at different hierarchical levels, where the raise in BB concentration during solvent evaporation displaced self-assembly equilibria toward more complex and differentiated structures. In addition, in a previous work, we determined that the presence of a free phenyl group at each of the two terminal ketimine moieties in structurally similar bis-salphen BBs made possible to activate their nanoscale molecular self-assembly into vesicles in solution, which upon solvent evaporation created 3D-networked structures that mimicked the neurons in a brain (Escárcega-Bobadilla et al., 2013a). In that work, those free phenyl groups gave access to further coalescence of the vesicles through drop casting to render interconnected microscale bodies at the mid-range mediated by π-π and van der Waals (vdW) interactions via variation of dielectric constant (ε) of the medium. In this work, we report on the design of a bis-Schiff-base molecular scaffold (e.g. molecule 1 in Figure 2) containing two key portions, α and β, which are connected through a central o-phenylenediamine derivative bridge, i.e. an unsymmetrical organic salphen scaffold. Portion α consisted of a similarly free phenyl group located at only one di-aromatic ketimine moiety, instead of two as reported earlier (Escárcega-Bobadilla et al., 2013a; Zelada-Guillén et al., 2018), while portion β represented a variable number of -OH groups at different relative positions in an aromatic aldimine moiety. Portion α was kept unchanged for all BBs (excepting for one of the blanks) to expect the role of an agonist component, which dictated positive feedback over mid-range aggregation of upper hierarchy level assemblies via π-π and vdW interactions regulated by ε of solvent. In this way, we pursued that the larger assemblies, due to their higher size and lower diffusivities in comparison with the free molecules, would favor local coalescence, thus playing as an agonistic signaling. Conversely, the β portion was varied to progressively control their access to intermolecular H-bonding depending on ε of aprotic solvents and, consequently, to promote an antagonic competition to the dispersive forces. In this manner, we expected a negative regulation of locally available matter for self-assembly through the highly diffusive free molecules able to migrate at longer ranges and, therefore, to render an antagonistic signaling.

Figure 2

Synthesized building blocks

In the common scaffold (e.g. 1–4), the di-aromatic ketimine moiety with a free phenyl group is kept unchanged at one part of the structure (part α, blue rectangle), while aromatic -OH groups at different positions in an aldimine moiety are incorporated at the other part of the molecule (part β, red circle); note that both parts are linked through an o-phenylenediamine derivative bridge. In blanks 5–6, the variable -OH at β is absent, while the phenyl group at α is present only for the former.

Results

Turing patterns by self-assembly of a single building block

Turing patterns (Figures 3A–3C) consisting of stationary waves with a constant λ centered at 1.45 μm, average range from 1.259 μm to 1.771 μm (Figure 3D), emerged at different levels of magnitude, when we drop cast acetone solutions of 1 at an onset concentration [1] = C (1 mM). The pattern was stably produced without relevant defects at surfaces that ranged from the micro (400 μm2 in Figure 3C) and the mesoscale (104 μm2 in Figure 3B) up to the macroscale (ca. 1 mm2 in Figure 3A but surfaces >1 cm2 were easily accessible), as observed by optical microscopy (OM). C was found to be the lowest concentration able to generate visible patterns after complete evaporation. Similar patterning results were also achievable for [1] up to 10 mM, but concentrations above and below this range resulted into irregular association at the microscale, i.e. random droplets or isolated microparticles, respectively. Turing patterns with different λ values were also obtained for acetone solutions from 2 and 3 at the same C concentration (Figures S1A and S1B). However, no pattern was obtained from neither 4 (Figure S1C) nor the blanks 5 (Figure S1D) and 6 (Figure S1E) in these conditions, which yielded only amorphous randomly distributed aggregates. Fast Fourier transform (FFT) of the respective OM data discarded any patterning from 4, 5, or 6, while FFT combined with linear grayscale analysis confirmed periodic waves from 1, 2, and 3 which consisted of isotropic populated zones (the crests) surrounded by unpopulated regions (the troughs) with an average distance λ between the crests that depended upon the -OH positions at β: cf. λ = 1.45 μm for 1 (Figures 3D and 3E), λ = 1.07 ± 0.327 μm for 2 (Figures S1F and S1K), and λ = 0.46 ± 0.081 μm for 3 (Figures S1G and S1l); no λ was measured for neither 4 (Figure S1H), 5 (Figure S1I) nor 6 (Figure S1J). Atomic force microscopy (AFM) confirmed the isotropic patterning trend for 1 (Figure 3F), where the analysis showed crest heights of ca. 60–100 nm (Figure 3G), wavelength confidence interval at p < 0.02 (CI) λ = 1.43 ± 0.172 μm, and crest widths between 0.8 μm and 1.5 μm; the data were consistent with FFT analysis in Figure 3G and grayscale analysis in Figure 3E. Confocal laser scanning microscopy (CLSM) showed that the crests in AFM and the populated waves in OM were internally composed of self-assembled vesicles (<1 μm diameter) that exhibited highly fluorescent spherical cores and non-emitting shells (Figure 3H), surrounded by remaining non-emitting amorphous matter. The vesicles followed a Gaussian distribution in a range between 0.4 μm and 1.05 μm, a mean value of 0.647 μm, and a standard deviation (s.d.) of 0.146 μm (Figure 3I). On the other hand, aprotic solvents with very different ε that could oppositely alter H-bonding capabilities in the molecules, such as chloroform or dimethyl sulfoxide (DMSO), produced either irregular droplets or amorphous microparticles for all the key structures when studied under comparable conditions (Figures 3J, 3K, and S1M–S1R). Hence, acetone resulted in a necessary input for the Turing patterning process to manifest during drop-casting experiments. Then, we used molecule 1 as a model system in this work to study such pattern formation from different approaches.

Figure 3

Turing patterns from a single molecule

(A) OM analysis of a drop-cast acetone solution [1] = C on glass shows the pattern at long-range in the macroscale; scale bar 100 μm, field: ca. 1 mm2.

(B) Analysis at higher amplifications shows the high stability of the pattern at the mesoscale; scale bar 10 μm, field: 104 μm2.

(C) In the microscale, OM shows a homogeneous distribution of dark waves (populated zones), separated by clear zones (unpopulated); scale bar 10 μm, field: 400 μm2.

(D) The 360° annular spot obtained by FFT analysis of B confirms the isotropic λ = 1.45 μm in the pattern (N = 171); range: 1.259 μm (N = 116) to 1.771 μm (N = 68).

(E) 1D grayscale analysis of the red line in C shows an example for the periodicity found at λ from FFT in D.

(F) AFM analysis of the specimens shows isotropic periodicity in crests.

(G) The height of the crests measured from the troughs in the red and cyan lines selected from F, raise up to 60–100 nm, show widths within the limit found in E (<1.5 μm) and λ = 1.43 ± 0.172 μm (confidence interval, CI, p < 0.02, ν = 23, under a t-distribution scheme).

(H) CLSM analysis shows that the crests are composed of self-assembled vesicles (<1 μm diameter) that are fluorescent at the inner part, which match the features found by OM and AFM; scale bar 20 μm.

(I) Particle size distribution analysis of vesicles in H yields a mean value of 0.647 μm following a Gaussian trend; standard deviation (s.d.) 0.146 μm, N = 147.

(J) OM micrograph from the molecule in chloroform at C shows irregularly spaced microparticles; scale bar 10 μm.

(K) OM micrograph from drop casting 1 at C in DMSO yields amorphous agglomerations; scale bar 10 μm. See also Figure S1 for OM from 2–6 at C = 1 mM in acetone (a–e), their FFT results (f–j), 1D grayscale analysis for 2 and 3 (k and l) and OM from 1–4 at C in CHCl3 and DMSO (m–r); scale bars in Figure S1 are 10 μm. Experimental details can be found in the STAR Methods section.

Identification of the morphogens and their roles in hierarchical self-assembly

We used dynamic light scattering (DLS), high-resolution transmission electron microscopy (HR-TEM), and 1H nuclear magnetic resonance (NMR) diffusion-ordered spectroscopy (DOSY), as well as data from molecular modeling to identify the components playing as the two morphogens A and B in acetone solutions of 1 and in their evaporated counterparts. First, DLS was used to study the self-assembly of 1 (Figure 4A) in solutions (j) at representative BB concentrations ranging from below to above the onset value C, so as to determine the population distribution of assemblies (k) and their hydrodynamic diameters (dh ) together with their diffusion coefficients (D) at different stages of the patterning process. Under steady-state conditions, the evolution in number of aggregates and their size as a function of [1] facilitated the identification of the population segment playing as the activator A. Solutions prepared at the onset concentration, j = C, and at 10-fold that concentration, j = 10C (a surrogate for a C sample after 90% of solvent evaporation), produced two different aggregate population segments, k = 1 and k = 2, according to our analysis of aggregate populations in solution; the first population component consisted of nanoscale assemblies (e.g. for C, dh  = 82 nm and D = 17 μm2 s−1), while the last one corresponded to microscale bodies of the same size order of the crest width observed by AFM and OM in dry counterparts (e.g. dh  = 1.63 μm). According to Turing’s mechanism during further stages in the patterning process, the difference in diffusion between both morphogens together with stochastic local perturbations which spontaneously caused that [A]local>>[B]local, also promote an adjacent regional depletion of A that results into a local transposition in the original concentration ratio between activator/inhibitor (i.e. from [A]0≥[B]0 to [A]local<<[B]local and [B]local<[B]0). This ratio inversion together with the local inhibitory effect driven by the predominance of B over A induces the more pronounced decrease in the relative local concentration of A in comparison with its initial steady-state value in the bulk before patterning occurred ([A]local<<[A]0), thus triggering formation of the void spaces and laterally propagating the instabilities until the final pattern is achieved. For a patterning system containing only one molecular component instead of two or more, where the oscillations produced in [A]local and [B]local must be mutually in-phase and should follow the Le Châtelier principle, such a reduction in [A]local is the microscale equivalent to diluting the activator in the macroscale. Therefore, we analyzed two different diluted solutions at [1]1] = 10−4 M). One of these solutions was prepared by diluting 10-fold the C stock solution (j = 0.1C) and the other one by dilution of the 10C sample in a factor of 1/100 (j = 10C), so that both reached equivalent concentrations under different routes. At these conditions, the microscale aggregates (k = 2) previously observed for j = C and j = 10C were now absent. Furthermore, the nanoscale assemblies (k = 1) reduced their average diameter in a factor of ca. 4 or 5, if respectively compared to their source solutions C or 10C. According to the model’s prerequisites, it is plausible that 1) at the onset conditions C, k = 1 nanoscale assemblies are playing the role of an active version of morphogen A, so that DA≡D = 17 μm2 s−1; 2) the k = 2 microscale aggregates coexisting at [1]≥C could be solvated components of the visible part in the pattern (a patterning readout prior to dryness), which might be produced by self-assembly of N units of morphogen A into a higher hierarchical level, following the Le Châtelier principle, the higher the concentration, the larger the self-association and vice versa; 3) from dimensional comparison of CLSM vs. DLS data, it is likely that k = 1 aggregates are the nanosized vesicles described earlier; 4) at diluted conditions [1]1] = 0.1C), the k = 1 nanoscale assemblies exist under a reversibly inhibited (inactive) form of morphogen A that does not trigger a patterning readout, which suggests that A should be progressively degraded during emergence of the void spaces. As free 1 molecules were neither confirmed nor discarded as a population segment by DLS, likely because of their low abundance, we used DOSY to confirm if they actually coexisted with the rest of the self-assembled counterparts (Figure S2). DOSY delivered a value of D = 1750 μm2 s−1 and dh = 8.0 Å, which according to density functional theory (DFT) and molecular dynamics (MD) calculations should correspond to the free molecule (Figure S3). Hence, despite the coexistence of morphogen A and solvated microscale pattern components at C, free molecules are the population component with the highest diffusivity in the pool (D>>DA, D ≈ 100DA) and thus, they are the most probable candidate for playing the role of the inhibitor morphogen B, so that DB≡D=1750 μm2 s−1.

Figure 4

Activation and inhibition of hierarchical self-assembly

DLS (A) shows morphogen A (level 2) in two reversibly stable forms in solution: inhibited A ([1]

DLS (A) shows morphogen A (level 2) in two reversibly stable forms in solution: inhibited A ([1]1]≥C) which triggers further self-assembly (level 3) into pattern components before dryness. HR-TEM confirms that A consists of discrete hollow spherical vesicle-type particles with a shell thickness τ (B) or solid structures (C) depending on the concentration; A can spontaneously degrade under dynamic conditions –i.e. bursting vs. reassembly in C–, through smaller nanostructured components (level 1); diameter of level 1 nanostructures (dL1) follows a Gaussian distribution (D), shaded orange area is ±1 standard deviation (s.d.), dL1 average is the arithmetic mean value, N = 48; scale bars are 20 nm for B and 50 nm for C. In silico studies (E) show that size d (spherical approximation), of n assembled molecules at level 1 evolves bi-logarithmically until a point where further self-assembly produce curved bodies (level 2) with thickness equivalent to τ in B; with the data, concerted activation/inhibition supramolecular reaction paths (F) are proposed from different hierarchies. The hydrodynamic radii of free molecule 1 at hierarchy level 0 were determined as 4.0 Å from DOSY using the Stokes-Einstein equation -Equation 1, this implies that the hydrodynamic diameter (dh) in solution is 8.0 Å, these data agree with the dimensions obtained as the minimum energy structure from the random generation of conformers followed by optimization at the DFT level, using the PBE(DTS)/DNP method with implicit solvation in acetone. Experimental and computational details can be found in the STAR Methods section. See also Figures S2–S4 and S11.

We used HR-TEM to gain insight on what should occur when the active morphogen A degrades to produce its inactive version (inhibition path). Drop-cast samples prepared at [1]<1] = 10−5 M) facilitated the disassembly processes and provided morphological information on the outcome of such path (Figures 4B–4D). However, as the inhibition path should be intrinsically connected to the activation path, in this section, we compared the experimental data from HR-TEM with the modeled size evolution of self-assemblies resulting from calculations (Figure 4E) to propose both paths under supramolecular reaction terms and hierarchical self-assembly (Figure 4F). The specimens analyzed by HR-TEM exhibited two types of spherical morphologies (cf. Figure 4B vs. Figure 4C) with diameter d values in agreement with DLS data for inhibited morphogen A (Figure 4A). The first morphology type consisted of hollow vesicle-type bodies (Figure 4B) with an external shell thickness, τ = 3.6–9.9 nm, where the smaller the value d, the thicker the shell (cf. Figure 4B vs. Figure S4). This trend was normally observed to a point where the structures became solid spherical nanoparticles for the smaller versions, while such tendency could on the contrary explain the existence of the larger fluorescent vesicles of ca. 600 nm diameter previously observed by CLSM in the final pattern at C (i.e. Figure 3H). Interestingly, we also found a prevalence of the solid spherical nanoparticles up to the highest concentration that in glass afforded the pattern ([1] = 10 mM, Figure S4). Occasionally, some of these inhibited morphogen A structures spontaneously carried out disassembly through a bursting process (Figure 4C). When the latter occurred, smaller nanostructures with a diameter of 4.1 ± 0.81 nm (Figure 4D) were expelled to the surroundings. The evolution of diameter as a function of the number of self-assembled molecules n (Figure 4E) evaluated through in silico studies indicated that these expelled nanostructures should be composed of ca. 8–21 molecules integrated in a hierarchical self-assembly level 1 (L1) in an inhibited version (namely L1). However, as shown in Figure 4C, it is likely that the disassembly of morphogen A into L1 nanostructures proceeded under dynamic conditions, in which L1 structures could reassemble again into larger spherical bodies of the same size order than morphogen A. Such scenario indicated that morphogen A should be structured at a higher hierarchy level 2 of self-assembly (i.e. L2). The latter was also supported by in silico studies, size evolution of modeled self-assemblies proceeded bi-logarithmically for oligomers between 2 and 256 molecules (Figure 4E), following an isotropic aggregation trend via surface minimization. In contrast, self-assembly of bodies larger than 256 molecules produced anisotropic and curved structures beyond this point. These curved assemblies were integrated by discrete numbers of L1 units in their activated version (namely L1), and their morphological features and thickness values matched the τ values measured by HR-TEM (cf. τ in Figure 4B vs. n = 1024 at L2 in Figure 4F).

These results strongly suggest that the external shell in morphogen A vesicles is composed of self-assembled L1 units if these reach oligomerization beyond a critical value. Consequently, the activation path for the patterning process should consist of both, oligomerization of free molecules at the hierarchy level 0 (namely L0) into L1 up to such a critical value (L1 nanostructures) and self-assembly of the latter into morphogen A at L2, through independent but concerted out-of-equilibrium supramolecular reactions successively favored by an increase in [BB] at a local level (Figure 4F, reaction path: a→a). Hence, as shown earlier by DLS data and according to Turing’s model, when the activation path occurs locally during solvent evaporation, the patterning process is triggered at the short- to the mid-range mediated by the assembly of morphogen A particles (the vesicles in Figure 3H) into solvated pattern components at a hierarchy level 3 (Figure 4A). On the contrary, when the inhibition path is favored, concerted dynamic equilibria should occur via degradation of L2 morphogen A into L1 nanostructures, until these oligomers finally may disassemble into free molecules at L0, i.e. in Figure 4F, reaction path: i→i. These L1 structures should actually be smaller than their L1 homologs involved in the opposed path (Figure 4F, products in i and a), such that their respective degrees of oligomerization cover size domains which do not overlap: cf. in Figure 4E, regression curve sections n:8 ≤ n ≤ 21 vs. n:32 ≤ n ≤ 256. In any case, as seen from HR-TEM Figure 4C, bursting of morphogen A into L1 nanoparticles might proceed dynamically to yield reassembled bodies composed of the latter, whereas infrequent totally scattered L1 structures seemed markedly reduced in size (apparently unstable). Then, both paths might be also connected through an intermediate equilibrium between the smaller oligomers L1 and the larger assemblies L1 (i.e. in Figure 4F, reactions i and a). In this way, depending on the local condition dynamics during the pattern emergence, different path combinations of these supramolecular reaction steps either at the activation path (a, a or a) or the inhibition path (i, i or i) could simultaneously occur in the system, thus creating a local equilibria shift that facilitates flux-balance, total mass conservation, non-linearity in reactions, and cooperativity during RD (Brauns et al., 2020, 2021; Diambra et al., 2015).

In silico studies in self-assembly modulated by supramolecular signaling

Molecular modeling considering the solvent ε using MD, DFT, molecular mechanics (MM), and Monte Carlo methods was carried out to identify the supramolecular mechanisms that rule the formation of self-assembled components in solution (supramolecular signaling). Through random binary mixing, we generated ten aggregate generations, up to 210 molecules each, finding that bodies composed of 1024 molecular units in acetone approximate the dimensions and morphology of the observed assemblies in morphogen A, as already mentioned in the last section. Plausible hierarchical levels L1 and L2 of self-assembly emerged during this expedition and are discussed concerning their formation.

In general, molecular modeling determined that depending on ε, H-bonding has paramount importance in self-assembly and seems to be decisive during the early growth of particles from free molecules in solution, as inferred from the three possible scenarios for self-assembly. This interaction may, however, be disrupted by highly polar solvents such as DMSO (ε = 46.7). In these conditions, vdW and π-π through portion α in the BB are the main interactions among molecules, with scarce H-bonding, as delivered from an analysis of the interactions in DMSO (Figure S5). The plasticity of these soft contacts, along with the likely occlusion of solvent, as well as the lack of further contribution to overall stabilization when particles grow beyond 32 molecules (Figures S6 and S7), are among the factors that can be associated to the absence of the Turing pattern using DMSO as solvent. On the other hand, H-bonding is especially important in low polarity medium, such as CHCl3, whose dielectric constant (ε = 5.62) exerts only a faint influence on the electrostatic component of H-bonding, as seen in our analysis of the interactions in chloroform (Figure S8). In these conditions, aggregation can be expected to be mainly led, at the first stages, by H-bond-driven oligomerization through portion β.

A moderately polar solvent such as acetone (ε = 20.7) proved to be key in tuning assembly of the bis-Schiff bases explored. Its polarity modulates H-bonding at the short-range, while its vapor pressure facilitates its departure during evaporation. At the molecular level, self-assembly of the most stable conformations of 1 occurs via partially hindered H-bonded domains through portion β, which appear as cooperative modes of oligomerization during earlier self-assembly generations in hierarchy L1 (Figures 5A and S9), with contribution from diffuse interactions at portion α such as π-π and vdW, as evidenced from an analysis of the interactions in acetone. This aggregation mechanism is less sterically hindered for Schiff bases 2 and 3, thus facilitating oligomerization through more accessible intermolecular H-bonding at portion β in these molecules (Figure S10). Therefore, accessibility of -OH groups at β might play a major role in the reduction of λ values for the patterns obtained from 1–3 (i.e. λ>λ>λ). The first level of organization L1 of 1 is likely H-bonding-driven assembly of small oligomers. During evaporation of relatively polar acetone, H-bonding is screened until concentration reaches a critical value, where the aggregation hierarchy L1 (Figure 5A) is a synergistic interplay of polar and dispersive contributions, modulated by ε and molecular steric hindrance. Higher order aggregates at L2 (Figure 5B) form through a combination of shape complementarity, diffuse vdW and π-π interactions from exposed α domains, as well as solvophobic interactions, spontaneously shifting from a homogeneous mixture of discrete smaller nanoparticles at L1 to the characteristic curved shell components of morphogen A at L2 (Figure S11). A similar process is expected to occur at least in part for 2 and 3; as assembly of aggregates proceeds (Figures S12 and S13), larger clusters in the order of ca. 103 molecules and up to 30 nm in size were formed under the modeling conditions tested, as shown during the analysis of the assembly dimensions in acetone. An estimation of aggregation and solvation energies on kinetically trapped supramolecular clusters (Figure 5C) determined that below a critical value starting at n = 8, solvation energy per molecule outweighs aggregation, and the dominant species in solution is L0 with infrequent occurrence of low self-assembly generation intermediaries (2 ≤ n < 8); this scenario should be favored under low local [1] conditions (e.g. during emergence of the void spaces in the pattern, when [B]local<[B]0 and [A]local<<[B]local). Above n = 32 (L1 onward), the aggregation energy per molecule has practically converged; this situation is expected to occur at higher local [1] regimes (e.g. during the formation of the visible parts of the pattern, when [A]local>>[B]local and [A]local>>[A]0). At L2 (n > 256), the solvation energy converges to values close to 0 kcal/mol; at this point, the tangent to each curve is practically zero and the assemblies begin behaving as pure phase. L1, L1, and L2 belong to the regime where aggregation energy overcomes solvation. L1 is closer to the threshold and is the transition zone where the change in particle size translates in a large change in solvation/aggregation energy differences (the slope of each curve is significantly ≠ 0 but aggregation is already exothermic). Therefore, volatility, which dissipates the overall exothermic contribution from self-assembly, and thus, could be seen as a type of energy currency, as well as polarity of the medium, are the two key factors to achieve total control of the desired assembly. These factors participate together with a regulated hindrance to intermolecular H-bonding through -OH groups at portion β of free molecules and during earlier generations of L1 self-assembly to rule the inhibitory factors in the patterning process via an antagonistic participation of H-bond signaling (information transfer) executed by highly diffusible components. On the contrary, the same two key factors also contribute to the activating component in the development of the patterns when self-assembly is driven by vdW and π-π at hierarchy L2 thanks to the further access to superficial aromatic domains provided by portion α of later generations at L1 through an agonistic involvement of vdW/π-π signaling carried out by little diffusible structures.

Figure 5

Supramolecular signaling at hierarchy levels

Molecular modeling of 1 in acetone determined that depending on the degree of oligomerization (n), self-assembly is signaled by either the portion β (at lower degrees) or α (at higher degrees). Earlier self-assembly generations (A) at L1 (e.g. L1) are mainly formed through regulated access to H bonding between the -OH groups at β by partially hindered H-bonding domains, according to particle modeling at the DFTB+ level with the 3ob Slater-Koster library; L2 (B) is an assembly of L1 units at later self-assembly generations (i.e. L1) via vdW/π-π interactions between domains containing exposed α portions. As self-assembly takes place (C), the weighed aggregation energy per molecule (♦) decreases steadily from 0 kcal/mol at L0, converging toward a nearly constant value of −13.2 kcal/mol per molecule for n ≥ 32 (L1 onward); solvation energies per molecule (▲) follow the opposite trend, they start from a minimum value for L0 (−20.6 kcal/mol) and consistently increase; intersection occurs at n ≈ 8; as a result, assembly is unstable for 2 ≤ n < 8. See also Figures S5–S13 and S22, Tables S2 and S4.

High-throughput mathematical analysis and the mechanism of Turing patterning

We used linear stability analysis to perform a high-throughput mathematical screening of the possible RD interaction networks (Figure 6A) between the different diffusible components at hierarchy levels L0, L1, and L2 (Marcon et al., 2016). This analysis was carried out to narrow down all the possible RD topologies associated with each network which are able to generate Turing patterns in terms of network feedbacks rather than reaction parameters, as well as to confirm if the roles proposed for A (L2) and B (L0) matched Turing’s model requisites. In this way, we expected to finally simulate what type of Turing patterns (spots, stripes, or labyrinths) would be expected from components with diffusivities as those measured so far, and if the analytically predicted patterning readout agreed with the experimental observations. In a first scenario, the mathematical screen was constrained with known diffusivities from experimental data for morphogens A (DA) and B (DB) for a 2-node RD network. In this manner, we could analytically assign the topology with the highest robustness, with the aim of finding the one with the highest volume of the kinetic parameter space that affords pattern-forming conditions with respect to all the universe of possible kinetic constant values. Such restriction was also done for a more complex 3-node RD network scenario also including the diffusivity approximations for the observed L1 hierarchy (DL1), using, respectively, the activated (L1) and inhibited (L1) versions as lower and upper interval limits for the range 107–346 μm2 s−1 (assumptions and data used for the estimations are indicated in the STAR Methods section). In this way, plausible supramolecular reaction mechanisms which are based only on the diffusible pool components L0, L1/L1, and L2 and the reaction equilibria from Figure 4F could be proposed for both 2-node and 3-node scenarios in order to define the RD role of each hierarchy and, at the same time, to determine the most probable mechanistic manifold in Turing patterning.

Figure 6

High-throughput mathematical analysis constrained to diffusion parameters

Automated analysis of all the possible 2-node and 3-node interaction networks (A) between hierarchy levels playing the morphogen roles A and B (black) or between A, B, and level 1 (black + gray), respectively delivered RD topologies B and C when constrained to the estimated diffusivities DA, DB, and DL1; topologies are consistent with experimental and modeling data, they confirm the roles for each hierarchy under an activator-inhibitor system and predict decay loops that are coherent with out-of-equilibrium supramolecular reactions in Figure 4F; κ are kinetic constant predictions (STAR Methods section, supramolecular reaction mechanisms for 2-node and 3-node RD networks). 1D simulations of B (D) and C (E) yield in-phase morphogen concentration spatial periodicities; 2D counterparts respectively predict Turing patterns F and G that are of the same type, similar dimensions, center for λ, and range values (see FFT insets) than the experimental results observed by OM and AFM in Figure 3; DL1 = 346 μm2s-1 for E/G calculations (DL1). See also Figures S14–S16. For Figure S16, 1D simulation (a) of the 3-node network produces morphogen A/B concentration periodicities which are in-phase with L1; 2D counterparts produce a Turing pattern (b) similar to experiments, but at a 7-fold larger wavelength; these results were obtained from DA = 17 μm2s-1, DB = 1750 μm2s-1, and DL1 = 107 μm2s-1. The 3-node RD network that resulted from predictions was identical to Figure 6C. Further details can be found in the STAR Methods section.

In both RD scenarios, topologies with reasonable positive and negative feedbacks resulted from the analysis (Figures 6B and 6C). These topologies were able to produce in-phase periodic variations in component concentrations (Figures 6D and 6E), while Turing patterns of the same family, dimensions, and λ intervals than the experimentally observed ones emerged after 2D numerical simulations (Figures 6F and 6G). In the first 2-node network scenario (up to four possible interactions, λ = 1.1 ± 0.44 μm), mathematical screening using the experimental data constrained RD topologies to a robust and physically coherent activator-inhibitor model candidate (Figure 6B) where a short-range activation was carried out by morphogen A, which was responsible for the periodic concentration crests in Figure 6D and the occupied red regions in the Turing pattern prediction Figure 6F. In this same scenario, a long-range inhibition was performed by B, where this morphogen was responsible for the periodic troughs in Figure 6D and the unoccupied black regions in the Turing pattern Figure 6F. This analysis also predicted a self-enhancing loop on node A that corresponded to its auto-catalyzed production as well as a self-regulation loop on B, which represented its self-clearance from the diffusible pool. Each feedback might be individually explained by proposed reaction mechanisms (Figure S14 and STAR Methods section, “supramolecular reaction mechanisms for 2-node RD network”) which involve reaction routes that include combinations from one to three of the steps a, a, a, i, and i in Figure 4F, depending on the case analyzed. Although reaction i is left out of this 2-node scenario, these possible mechanisms are coherent with experimental data and molecular modeling results, whereas the predicted kinetic constant intervals estimated for each interaction could give place to robust 2D Turing patterns in the range for λ from 0.66 μm to 1.54 μm (1.45 μm experimentally). However, such topology is not able to provide a reasonable explanation for the diffusion role of L1 nanostructures on the overall patterning process. On the contrary, participation of L1 in the proposed mechanism reaction steps is restricted to an intermediary level despite the evidence of their role played in hierarchical self-assembly, thus limiting the scope of the full 2-node RD mechanism to a simplified scenario.

In this sense, mathematical analysis through a 3-node RD network scenario (up to nine possible interactions, λ = 1.8 ± 0.58 μm when L1 data is tested) offered a full picture and a more comprehensive mechanistic manifold (STAR Methods section, “supramolecular reaction mechanisms for 3-node RD network” and Figure S15) to the overall patterning process observed. Under the 3-node scenario, the resulting RD topology consists of a short- to mid-range activation process simultaneously carried out by the little diffusible pair morphogen A/hierarchy level 1 (L2/L1), where the activation information (the agonistic supramolecular signaling via vdW and π-π through portion α in 1) is hierarchically transferred from L2 to L1 (and vice versa) to finally deliver the same from L1 to morphogen B (positive feedbacks in Figure 6C). This process is concerted with a long-range inhibition driven by the highly diffusible morphogen B, which directly interacts with both activating nodes (negative feedbacks in Figure 6C), primarily via antagonistic H-bonding signaling through portion β in 1. However, L1 also exhibits a self-mediated degradation (negative feedback loop in Figure 6C) at the mid-range and thus, a self-imposed antagonistic signaling through similar H-bonding interactions as in B. In this topology, the activation and inhibition processes are also possible under both limiting scenarios for DL1, but L1 affords wavelength values 1.22 μm to 2.38 μm that are closer to the experiments than those from L1. Either L1 or L1, together with L2 morphogen A, synergistically impose the spatially periodic concentration crests in Figure 6E (for L1) and Figures S16A (for L1) as well as the occupied regions in the predicted Turing patterns (red regions in Figures 6G and S16B). Under both DL1 limit scenarios that respectively yield periodic wavelengths at λ = 1.8 ± 0.58 μm for L1 and λ = 10.7 ± 3.16 μm for L1, the concentration troughs in Figures 6E and S16A are ruled by a direct inhibitory effect of morphogen B on both upper hierarchy levels and the self-clearance predicted for L1. These factors are responsible for the unoccupied regions in the Turing pattern simulations (black regions in Figures 6G and S16B). The analysis also predicts kinetic constant ranges for each interaction which are more reasonably explained in terms of all the 6 reactions a and i (Figure 4F) through routes that involve from one to three steps. These steps could be classified into five different types of matter transfer processes which depend upon the reaction step in the mechanism: (self) assembly (a and a), coalescence (a), splitting (i and i), leaching, or sequestering (i, i, and i for both). However, these five matter transfer scenarios also depend upon the feedback class, where the sum of positive feedbacks excludes the participation of sequestering processes, while the negative feedbacks do not include (self)assembly. Evidence that suggests their viability in the overall proposed mechanism could be found in the results from experiments/molecular modeling/predicted spatial concentration oscillations. Moreover, this 3-node-based mechanism is more solidly sustained by the solvation vs. aggregation energy profiles if compared with the simpler 2-node mechanism.

Discussion

Our findings indicate that Turing pattern emergence is carried out with a synergistic and cooperative participation of the three diffusible nodes L0, L1, and L2, which trigger the patterning readout at the hierarchy L3. The supramolecular reaction mechanisms proposed for the positive feedbacks are meeting both, the Turing’s model, and predicted patterning dimensions according to the estimated diffusion values and proposed morphogen roles, predominantly ruled by agonistic supramolecular signaling based on vdW and π-π interactions mediated by little diffusible particles. On the contrary, the mechanism counterparts that explain the negative feedbacks are mediated by antagonistic signaling through H-bonding from highly diffusible pool components. It is therefore reasonable to expect that also the rate-determining step in each of the mechanistic routes is dominated by these supramolecular interactions depending upon whether they are part of a positive or a negative feedback. However, such a level of more complex understanding is at this point an open question and currently ongoing studies aim to unravel the roles of variable self-assembly rates arisen from the dynamics in size-evolution and hierarchical regulation of non-covalent interactions.

Nevertheless, it is important to highlight that the spontaneous formation of Turing patterns by hierarchical self-assembly of a single molecular building block has never been reported so far. Therefore, this study shows that it is possible to use differentially modulated supramolecular interactions as either agonistic or antagonistic signaling regulators in order to transfer molecular chemical information up to the mesoscale in the form of a patterning readout. This situation opens evidence-based conceptual venues toward more complex patterning mechanisms departing from only one elementary assembling unit, where the classical yet simple 2-node activator-inhibitor RD system has practically remained the basis to explain RD Turing patterns from different reactants for six decades. In this scenario, it is now much easier to hypothesize on, for example, how life evolved from inert molecular components using not only chemical transformations but also supramolecular interactions as a currency for information exchange, which in turn, could more efficiently be exploited by protocells in different adaptive manners depending upon the number of entities associated in a constrained space.

Limitations of the study

In this study, we only identified the minimum molecular programming requisites for an organic salphen compound to drive self-assembly toward simultaneous populations of different diffusivities and sizes in acetone solutions, which upon evaporation rendered a Turing pattern in the solid state. We explained the process hierarchically, through concerted and networked supramolecular reactions based on the solely constructive or destructive matter interchange scenarios allowed for the case. However, we identified two plausible RD topologies that involved activation and inhibition expressed in the form of positive and negative feedbacks under kinetic ranges that were only predicted analytically. We circumvented this limitation by deploying not only a comprehensive molecular dynamics study coupled to the mathematical analysis but also achieved the largest number of experimental data possible for such a dynamic heterogeneous system, from OM, CLSM, AFM, HR-TEM, DLS, and NMR DOSY.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, Gustavo A. Zelada-Guillén (g.zelada@unam.mx).

Materials availability

All unique reagents generated in this study are available from the lead contact without restriction.

Experimental model and subject details

This study did not use experimental models typical in life sciences.

Method details

General synthesis and characterisation

All starting materials were purchased from commercial sources and used without further purification. Compounds (E)-2-(((2-aminophenyl)imino) (phenyl)methyl) phenol (Escárcega-Bobadilla et al., 2012) and the respective blank BBs 2-((E)-((2-(((E)-2-hydroxybenzylidene)amino)phenyl)imino) (phenyl)methyl)phenol, 5 (Escárcega-Bobadilla et al., 2013b), and 2,2’-((1E,1′E)-(1,2-phenylenebis(azaneylylidene))bis(methaneylylidene))diphenol, 6 (Deng et al., 2014), were prepared using reported methodologies as described ahead; briefly, stoichiometric amounts (ca. 0.35 mmol each, ratios 1:1 for first and second compounds, 2:1 for the third one) of the necessary carbonylic synthon (e.g. 2-hidroxybenzophenone for the first compound, 2-hydroxybenzaldehyde for the second and third ones) and the amine precursor (o-phenylenediamine for the first and third compounds, (E)-2-(((2-aminophenyl)imino) (phenyl)methyl) phenol for the second one) were dissolved in 10 mL of methanol (MeOH) and stirred for 18 h at 298 K; the solids were filtered and washed with MeOH and their identities were confirmed by NMR using a suitable deuterated solvent as described ahead. Elemental analyses (EA) were performed with a Perkin Elmer 2400 for CHNS analyser. All NMR measurements were carried out on a Varian VNMRS 400 MHz apparatus at 298 K in acetone-d6 or CDCl3 as needed, and chemical shifts are given in ppm vs. TMS. High resolution mass spectrometric (HR-MS) data were obtained with an Agilent 6530 QTOF spectrometer from Unidad de Servicios para la Industria Petrolera (USIP) of the School of Chemistry – UNAM. DOSY experiments were carried out on a Bruker Ascend 500 MHz apparatus at 298 K in acetone-d6 and [1] = C and were determined at Laboratorio Universitario de Resonancia Magnética Nuclear (LURMN) – UNAM. The hydrodynamic radii of 1 was determined as 4.0 Å from DOSY results using the Stokes-Einstein equation -Equation 1-, this implies that the hydrodynamic diameter (dh) in solution is 8.0 Å; see also Figure S2.

Synthesis of bis-schiff-base 1

3-((E)-((2-(((E)-(2-hydroxyphenyl) (phenyl)methylene)amino)phenyl)imino)methyl)benzene-1,2-diol, 1, was synthesised as described ahead. (E)-2-(((2-aminophenyl)imino) (phenyl)methyl) phenol (100 mg, 0.35 mmol) and 2,3-dihydroxybenzaldehyde (51 mg, 0.37 mmol) were dissolved in 10 mL of MeOH and refluxed at 337 K for 18 h. A red suspension was formed, the solid was filtered and washed with MeOH. Red solid, 140 mg, 98%. 1H NMR (400 MHz, Acetone-d6): δ = 6.79 (ddd, JHH = 8.2, 7.4, 1.2 Hz, 1H), 6.83 (t, JHH = 7.8 Hz, 1H), 6.95–6.98 (m, 2H), 7.01–7.05 (m, 2H), 7.06 (dd, JHH = 7.8, 1.6 Hz, 1H), 7.1–7.16 (m, 2H), 7.18–7.20 (m, 2H), 7.27–7.36 (m, 4H), 7.43 (ddd, JHH = 7.8, 1.6 Hz, 1H), 8.67 (s, 1H, CHaldimine), 8.67 (s, 1H, OH), 13.25 (s, 1H, OH), 13.94 (s, 1H, OH). 13C{1H} NMR (100 MHz, CDCl3): δ = 175.11, 162.71, 161.94, 149.31, 145.28, 142.00, 139.25, 134.65, 133.69, 132.55, 129.12, 128.48, 128.04, 127.33, 125.69, 123.34, 122.94, 119.72, 119.11, 118.85, 118.22, 117.80. MS (MALDI +) m/z 408.732 ([M+]). EA calcd. (%) for C26H20N2O3 C 76.46, H 4.94, N 6.86 found C 76.23, H 4.66, N 6.81. For 1H and 13C{1H} NMR spectra see Figure S17.

Synthesis of bis-schiff-base 2

4-((E)-((2-(((E)-(2-hydroxyphenyl) (phenyl)methylene)amino)phenyl)imino)methyl)benzene-1,3-diol, 2, was synthesised as described ahead. (E)-2-(((2-aminophenyl)imino) (phenyl)methyl) phenol (100 mg, 0.35 mmol) and 2,4-dihydroxybenzaldehyde (50 mg, 36 mmol) were dissolved in 10 mL of MeOH and refluxed at 337 K for 18 h. A red suspension was formed, the solid was filtered and washed with MeOH. Red solid, 136 mg, 95%. 1H NMR (500 MHz, Acetone-d6): δ = 6.32 (d, JHH = 2.3 Hz, 1H), 6.48 (dd, JHH = 8.5, 2.3 Hz, 1H), 6.79 (td, JHH = 7.6, 7.2, 1.2 Hz, 1H), 6.95 (dt, JHH = 5.5, 3.8 Hz, 1H) 7.03 (ddd, JHH = 7.6, 6.0 Hz, 1.5), 7.10 (dd, JHH = 6.0, 3.4 Hz, 1H), 7.21 (d, JHH = 6.6 Hz, 2H), 7.26–7.36 (m, 4H), 7.40 (d, JHH = 8.5 Hz, 1H), 7.43 (ddd, JHH = 8.5, 7.2, 1.5 Hz, 1H), 8.61 (s, 1H, CHaldimine), 9.14 (s, 1H, OH), 13.49 (s, 1H, OH), 14.02 (s, 1H, OH). 13C{1H} NMR (100 MHz, CDCl3): δ = 175.76, 164.60, 163.50, 163.27, 163.11, 142.83, 135.36, 134.24, 133.05, 129.79, 129.08, 128.95, 128.76, 127.18, 126.42, 123.65, 120.55, 118.87, 118.79, 118.61, 117.28, 113.76, 108.53, 103.43. MS (MALDI-TOF) m/z found 389.91 ([M+−H2O]). EA calcd. (%) for C26H20N2O3·¼MeOH C 75.71, H 5.08, N 6.73 found C 75.30, H 4.89, N 6.72. For 1H and 13C{1H} NMR spectra see Figure S18.

Synthesis of bis-schiff-base 3

4-((E)-((2-(((E)-(2-hydroxyphenyl) (phenyl)methylene)amino)phenyl)imino)methyl)benzene-1,2,3-triol, 3, was synthesised as described ahead. (E)-2-(((2-aminophenyl)imino) (phenyl)methyl) phenol (100 mg, 0.35 mmol) and 2,3,4-trihydroxybenzaldehyde (57 mg, 35 mmol) were dissolved in 10 mL of MeOH and refluxed at 337 K for 18 h. A brown suspension was formed, the solid was filtered and washed with MeOH. Brown solid, 146 mg, 98%. 1H NMR (500 MHz, Acetone-d6): δ = 6.49 (d, JHH = 8.5 Hz, 1H), 6.79 (ddd, JHH = 8.3, 7.3, 1.2 Hz, 1H), 6.92–6.93 (m, 1H), 6.97 (d, JHH = 8.5 Hz, 1H), 7.02 (dd, JHH = 8, 1.6 Hz, 1H), 7.03 (dd, JHH = 8.3, 1.2 Hz, 1H), 7.07–7.11 (m, 2H), 7.18–7.20 (m, 2H), 7.24–7.25 (m, 1H), 7.28–7.35 (m, 3H), 7.43 (ddd, JHH = 8.6, 7.2, 1.7 Hz, 1H), 7.72 (s, 1H, OH), 8.24 (s, 1H, OH), 8.53 (s, 1H, CHaldimine), 13.48 (s, 1H, OH), 14.01 (s, 1H, OH). 13C{1H} NMR (101 MHz, Acetone-d6): δ = 175.93, 163.85, 163.53, 152.03, 150.83, 142.61, 140.78, 135.53, 134.30, 133.15, 133.03, 129.86, 129.04, 128.81, 127.20, 126.54, 125.35, 123.75, 120.62, 119.25, 118.87, 118.68, 113.94, 108.47. MS (MALDI+) m/z 425.462 ([M+H]). EA calcd. (%) for C26H20N2O4 C 73.57, H 4.75, N 6.60 found C 73.31, H 4.51, N 6.72. For 1H and 13C{1H} NMR spectra see Figure S19.

Synthesis of bis-schiff-base 4

2-((E)-((2-(((E)-(2-hydroxyphenyl) (phenyl)methylene)amino)phenyl)imino)methyl)benzene-1,4-diol, 4, was synthesised as described ahead. (E)-2-(((2-aminophenyl)imino) (phenyl)methyl) phenol (100 mg, 0.35 mmol) and 2,5-dihydroxybenzaldehyde (50 mg, 36 mmol) were dissolved in 10 mL of MeOH and refluxed at 337 K for 18 h. A yellow suspension was formed, the solid was filtered and washed with MeOH. Yellow solid, 143 mg, 99%.1H NMR (400 MHz, Acetone-d6): δ = 6.73 (d, JHH = 8.8 Hz, 1H), 6.77 (ddd, JHH = 8.3, 7.2, 1.3 Hz, 1H), 6.92 (dd, JHH = 8.8, 3.0 Hz, 1H), 6.95–7.03 (m, 4H), 7.12 (ddd, JHH = 7.1, 4.7, 1.9, 2H), 7.19 (dt, JHH = 6.6, 1.7 Hz, 2H), 7.27–7.35 (m, 4H), 7.41 (ddd, JHH = 8.6, 7.2, 1.7 Hz, 1H), 8.02 (s, 1H, OH), 8.62 (s, 1H, CHaldimine), 12.43 (s, 1H, OH), 13.92 (s, 1H, OH). 13C{1H} NMR (101 MHz, Acetone-d6): δ = 174.77, 162.72, 161.98, 155.26, 148.05, 141.29, 139.69, 134.31, 133.65, 132.42, 129.05, 128.38, 128.20, 127.95, 127.06, 125.70, 123.20, 121.63, 119.44, 118.90, 118.42, 118.18, 118.11, 118.06, 117.33. HRMS (MALDI +) m/z 409.18981 ([M+H]). EA calcd. (%) for C26H20N2O4·½H2O C 74.81, H 5.07, N 6.71 found C 74.46, H 5.07, N 6.71. For 1H and 13C{1H} NMR spectra see Figure S20.

Optical microscopy

OM analysis was carried out with a Motic model BA210 microscope equipped with Live Imaging Module Digital Camera. For greyscale and FFT analysis, the images were processed and statistically analysed with ImageJ 1.50b software (Wayne Rasband). Samples were prepared by drop-casting 1–10 mM solutions (concentration deviations are mentioned in the text) of the corresponding compound (previously dissolved in the selected solvent) on a glass slide at 298 K and dried completely for 1 day. Fields were adjusted to highlight the representative features of the samples analysed. See Figures 3 and S1.

Confocal Laser Scanning Microscopy

CLSM analysis was carried out with an Olympus model FV1000 microscope, using a green diode laser excitation wavelength of 559 nm. The images were processed with FV10-ASW Ver. 01.07c software, (Olympus Life Science) and statistical analysis was processed with ImageJ from data collected during emission mapping. Samples were prepared by drop-casting 1 mM solutions of 1 (acetone) on a glass slide at 298 K and dried completely for 1 day. See Figure 3.

Atomic force microscopy

AFM analysis was performed at tapping mode on a Park Systems model NX10 microscope using an AC160TS cantilever and the data was processed with XEI package (Park Systems). Samples were prepared by drop-casting solutions of 1 (1 mM in acetone) on a glass slide at 298 K and dried completely for 1 day. See Figure 3.

Dynamic light scattering

DLS analysis was carried out on a Zetasizer Nano ZSP (Malvern Instruments Ltd.) with 532 nm laser radiation source. 3 mL of the corresponding solutions of 1 (variable concentration in acetone) were analysed on a standard fluorescence glass cuvette with an optical path length of 1 cm at 298 K. See Figure 4.

HR-TEM

High Resolution Transmission Electron Microscopy (HR-TEM) analysis was performed on a JEOL model JEM-2010 microscope with an acceleration voltage of 200 keV, 50 micrometer C2 aperture, spot size 1 and variable dose rate. Images were processed with Gatan Digital Micrograph software (Gatan, Inc.) and analysed with ImageJ. Liquid samples were prepared by drop-casting 10 μM to 10 mM solutions of 1 in acetone on formvar carbon film-covered square mesh copper grids and dried completely for 4 h at 298 K. See Figure S4. In the same figure, TEM analysis of drop-cast solutions from 1 in acetone (10−5 M upper micrographs, 10 mM central-inferior micrograph) shows that in general terms for the diluted regimes, the larger the diameter d in hollow spherical vesicle-type bodies (morphogen A), the shorter the thickness τ of the external shell, but solid nanosized spherical counterparts are prevalent even at the largest concentration tested; scale bar is 50 nm for the three micrographs in the upper row; scale bar is 400 nm for the central-inferior panel, field: 16 μm2. See also Figure 4.

In silico studies

Calculations were carried out using the Materials Studio 8 suite (MS8), provided by Biovia, Dassault Systèmes, unless stated otherwise. The DMol3 program was employed for DFT computations with the COSMO implicit solvation model as implemented in MS8. Solvation energy benchmarking was performed using the SMD method (Marenich et al., 2009) with the ORCA 4.0 program (Neese, 2012). See Figures 4, 5, S3, and S5–S13.

Solvation and aggregation energies on supramolecular clusters were semi-quantitatively estimated upon calibration of an MM-based methodology and their detailed heuristics are explained in further subheadings. Briefly, we proposed random binary mixing to approach the modelling of kinetically trapped conglomerates and explore the energetic aspects of the first steps of self-assembly leading to nanoparticle consolidation. To calibrate our MM method, intramolecular H-bonding and rotation of the -OH group were the first parameters to be benchmarked. Our reference values came from DFT calculations using the PBE functional (Perdew et al., 1996) with the Tkatchenko-Scheffler empirical dispersion scheme for geometry optimisations (Tkatchenko and Scheffler, 2009), and the double-zeta polarised basis DNP, PBE(DTS)/DNP, with final energies evaluated with the m06L functional and the triple-zeta TNP basis, in vacuo for direct comparison with MM. Using the ESP-fitted DFT-derived atomic charges and the Dreiding (Mayo et al., 1990) force field (DFF) we succeeded to locate a maximum along the dihedral scan and the endothermic character of the local minimum at 180° was qualitatively correct. This protocol was subsequently used for all tasks performed at the MM level.

A set of conformers was generated for each molecule with the following methodology: a diluted (1 molecule in 1000 units of solvent) solution in acetone, with density 0.78 g/mL, was built using the Amorphous Cell program and the MM protocol; systems were then subject to molecular dynamics (carried out with the Forcite code, covering 1000 ps in 1 fs steps, under the NVT ensemble at 473 K to facilitate conformational flexibility). From the resulting trajectories, five frames were extracted, at 200 ps sampling intervals. The molecules from each frame were stripped from solvent molecules and optimised.

Supramolecular assemblies were generated using the Blends program, an implementation of the molecular silverware methodology (Blanco, 1991). Stochastic binary mixing was modelled, generating random aggregates of 2x molecules, with x ranging from 0 to 10, meaning x the generation number in self-assembly so that 2 = n. Using this strategy, particles reaching 200–300 Å in length and up to 1024 molecules (52,224 atoms per aggregate) were built and closely inspected, studying their self-assembly from single molecules to nanoparticles. With this static approach, we maintained atomic resolution throughout the study while accessing particle sizes of up to 1,024 (210) molecules with reasonable computational cost.

Geometry benchmark of aggregates was then carried in two steps. (1) Optimisation of bimolecular mixtures using (a) the DFT-based tight binding method DFTB+ (Aradi et al., 2007) with the ob3 Slater-Koster library (Gaus et al., 2013) and (b) the DFT PBE(DTS)/DNP procedure. Structural similarity between these sets was quantified as the overlap of steric fields, defined as the Lennard-Jones potential seen by a probe carbon atom on each structure. The structural match between DFT and DFTB+ was on average 99 ± 1%, so this method proved robust for the generation of DFT quality geometries for a fraction of the cost. On step (2) structures optimised with the MM protocol were compared with those obtained through DFTB + for generation 1 (89 ± 5% similarity), generation 2 (85 ± 3%) and generation 3 (78 ± 2%). These deviations reflect the flat potential energy surface associated with larger aggregates yet are acceptable given the large scale of the calculations performed and, thus, we applied this MM protocol for the general geometry optimisation of molecules and clusters.

The aggregation energy benchmark was performed using a set of conformers which were generated using the Conformers program with the MM protocol. From the results, the lowest energy conformer was optimised and used as seed to generate a collection of 105 bimolecular mixtures with the Blends program. The 100 lowest energy pairs were kept for the next step. Without further optimisation, aggregation energy was computed at the DFT level using the m06L/TNP method in vacuo. By comparing average aggregation energies from both methods (DFT: −14.4 ± 1.2 kcal/mol, MM protocol: −13.5 kcal/mol), we conclude that DFF embodies an adequate parameterisation yielding semi-quantitative aggregation energies vs. DFT computed values.

Solvation energy calculations were benchmarked as follows. For a given geometry, solvation was computed at different levels: a) Reference: m062×/def2tzvp with the SMD method, calculated in ORCA 4.0; b) m06L/TNP with COSMO solvation in acetone; c) random mixing with the chosen force field and the Blends program in MS8. For single molecules (generation 0), COSMO was +1.5 kcal/mol off and the Blends method reproduced the SMD solvation energy within 1 kcal/mol. For generation 1 aggregates, COSMO was ca +8.3 kcal/mol off and Blends was +0.5 kcal/mol off. This random mixing approach yielded semi-quantitatively solvation energies, while being relatively inexpensive. Therefore, the Blends method was standardised to the following procedure: Perform three mixing runs. The average base + screen interaction energy EBS is stored; three coordination number runs with 30, 50 and 150 iterations are run to extrapolate coordination numbers, C. Solvation energy is obtained as ΔEsolv = EBS x C. The method proved affordable up to generation 7 (n = 128). Solvation for larger aggregates was obtained through extrapolation leveraging on the linear relationship between size and solvation energy found for generations n > 4.

High-throughput mathematical analysis

Automated linear stability analysis and high-throughput analytical prediction of Turing patterning conditions were carried out using the web-based software RDNets (Marcon et al., 2016) through the Wolfram CDF Player 11.3 interface (Wolfram Research). The algorithm used has been recently reported to successfully predict Turing pattern-forming diffusivity and kinetic constant parameter ranges in RD networks by exploiting a new graph-theoretical formalism that simplified analysis to a graphical-user interface under a reasonable computational demand. See Figures 6 and S16.

During the mathematical analysis, we left unconstrained all the possible interactions (activation or inhibition feedbacks) between each hierarchy, the pattern phase (in-phase, out-of-phase) and the pattern type (I, II or III) in order to avoid bias in the final results. As a first step, we constrained the 2-node RD scenario with the two morphogen diffusivity values DA = 17 μm2s-1 for L2 and DB = 1750 μm2s-1 for L0; for the 3-node scenario, we used the same previous hierarchy diffusivities together with a third node representing L1, in which we alternatively included either the lower limiting value DL1 = 107 μm2s-1 for L1 or the upper limiting value DL1 = 346 μm2s-1 for L1, to consider a dynamic case in which L1 can progressively cover a range of sizes rather than comprising size-defined independent nodes, and still be able to participate in the Turing patterning process.

Afterwards, from the different RD networks generated through each mathematical analysis, we only selected the systems with the highest estimated robustness per interaction network and discarded the rest. These selections were further analysed through the software pipeline to determine their respective network topologies in order to predict the kinetic constant intervals (κ>0 for the activation feedbacks, κ<0 for the inhibition feedbacks, 2-node:1 ≤ l ≤ 4, 3-node: 5 ≤ l ≤ 10) that delivered viability in Turing patterning under the diffusivity constrains provided at the beginning of the analysis; for computational reasons, the 3-node scenario was analysed in two steps: first, prediction of κ to κ general conditions was carried out automatically, and second, prediction of κ was performed using selected data from the ranges given by the former step (selected data is described ahead). Subsequently, by using representative numerical data for each κ (×10−4 s−1) and variable final times (s), we carried out 1D analysis of the concentration fluctuations for the diffusible components vs. space: for 2-node analysis we used κ = 0.504, κ = 1, κ = 1, κ = 1, final time 61.951; for 3-node analysis with L1 we used κ = 1, κ = 1, κ = 0.00391, κ = 1, κ = 0.00171, κ = 0.00915249, final time 15.530; for 3-node analysis with L1 we used κ = 1, κ = 1, κ = 1, κ = 1, κ = 0.625, κ = 0.1875, final time 124.24. Finally, we used the same feedback constrains to predict spatial 2D Turing patterning capabilities for morphogen A concentration under each scenario, to determine whether the graphical results were consistent with the labyrinth patterning families observed in the experiments.

Analysis of aggregate populations in solution

For j = C solutions, k = 1 consisted of nanoscale assemblies with dh  = 82 nm (D = 17 μm2s-1, from the Stokes-Einstein equation -Equation 1-), whereas k = 2 corresponded to a 20-fold larger aggregate at the microscale with dh  = 1.63 μm (D = 0.86 μm2s-1), with a volume ratio k = 2:k = 1 of V ≈ 7.8×103V under a spherical model. Interestingly, the microscale population presented the same size order of the crest width observed by AFM and OM in dry counterparts.

We also analysed a ten-fold concentrated solution (j = 10C), as a surrogate for a C sample after 90% of solvent evaporation, in order to determine what occurs with the aggregate populations during the drop-casting process but before dryness. The latter produced a simultaneous increase of dh in a factor of ca. 3 for both populations k = 1 and k = 2, giving larger nanoscale and microscale counterparts with dh  = 227 nm and dh  = 4.56 μm. In the 10C solutions, coexistence of both aggregates k = 1 and k = 2 (D = 6.2 μm2s-1 and D = 0.31 μm2s-1, respectively) occurred as in the original C solution. These populations also preserved the same diameter ratio 1:20 and roughly the same volume ratio (V ≈ 8.1 × 103V) as in C, but respectively represented an increase in spherical volume equivalent to 21–22 times that resulted from rising the concentration ×10 (cf. V ≈ 21V vs. V ≈ 22V).

In line with the observed trend, if k = 2 is produced from N self-assembled k = 1 particles through a non-crystalline spherical agglomeration model with a lack of inter-particle spacing due to plasticity in the components, then their volume ratio would be a fair approximation to the number of k = 1 units integrated in one k = 2 aggregate. Therefore, in either C or 10C, microscale bodies k = 2 are derived from the coalescence of ca. 7.8 × 103–8×103 nanoscale particles (N) of the population k = 1 observed at each concentration.

On the other hand, their diluted counterparts j = 0.1C and j = 10C did not favour the expression of microscale populations k = 2, but only produced nanoscale assemblies k = 1 with 4-fold or 5-fold decrease in diameter values: dh  = 21 nm (D = 67 μm2s-1) or dh  = 44 nm (D = 32 μm2s-1).

Three possible scenarios for self-assembly

Self-assembly in the studied bis-Schiff bases can be divided in three scenarios, depending on the polarity of the solvent. In the first scenario, a polar solvent with a high H-bond disrupting character, such as DMSO, rules out oligomerisation through H-bonding; hence, short range contacts such as π-π and vdW interactions are expected to dominate self-assembly. In the second one, involving a low-dielectric constant medium such as CHCl3, the electrostatic component of H-bonding is key to aggregation. This interaction is not screened and thus may lead to the anisotropic formation of H-bonded oligomers. In the third case, acetone, being intermediate between these extreme cases, can be expected to present characteristics from both the first and second scenarios.

In the computational study, factors such as the evaporation time of solvents (1.9 for CHCl3, 1.8 for acetone and 1500 for DMSO, considering Et2O evaporation as unity) and their viscosity (0.57, 0.33 and 2.0 cP, respectively) -Smallwood, 1996-, are not accounted for, yet these factors are expected to add to the full picture of the complex phenomenon under analysis.

Analysis of the interactions in DMSO

In highly screening DMSO, H-bonding can be regarded as nearly absent. This leads to the most isotropic scenario, where no directional forces influence aggregation. The short-range of π-π and vdW interactions means that, in these conditions, coalescence is a random process where a distribution of configurations interact upon contact without preferential modes. We simulated this by random binary mixing, where stable conformations are blended without bias.

A close look to the sequential growing of aggregates allows delving further into the manifold of plausible self-assemblies. Ideally, first-generation aggregates would be bimolecular with suppressed influence of intermolecular H-bonding. Due to the presence of 4 aromatic rings, dispersive forces are the main contributions to intermolecular contacts, especially parallel-displaced and edge-to-face π-π and vdW interactions. As estimated using the DFT methodology, bimolecular interactions are energetically close-lying if H-bonding is excluded. As clusters grow, the gamut of intermolecular contacts broadens and the appearance of accidental H-bonding is expected. As stated by Day and co-workers in their studies on molecular crystals (Thompson and Day, 2014) and herein applied to nascent nanoparticles, the contacts which maximise surface overlap can be strongly favoured even at the cost of conformational strain. A rough visual analysis of Figure S5 shows that the most compact aggregates are those where the interacting surface of the conforming molecules are maximal. This can be illustrated with the surface of the aggregate that is accessible to solvent. As will be shown, this approximate relationship holds true from the bimolecular cluster space up to particles in the nanoscale, composed of > 103 molecules. See also Figure 5.

In Figure S5, self-assembled structures are obtained when the influence of intermolecular H-bonding in 1 is partially suppressed by a high dielectric constant scenario and no solvent occlusion occurs. First-generation aggregates are bimolecular, second and third generations are particles constituted by 4 and 8 molecular units, respectively. Aggregates include aggregation energies and solvent-accessible surfaces for each particle. Energy values shown are relative to the free units and weighted with the number of molecules in the cluster, reflecting the average aggregation energy per molecule.

Particle surface naturally increases with the number of constituting units; this value, weighted with the number of molecules in the aggregate, shows the trend observed in Figure S6. As aggregates grow near 32 molecules, the surface/size ratio becomes nearly constant; further coalescence does not contribute to decrease this ratio. If the total and weighted aggregation energies, rather than surface, are now compared in a similar fashion, it can be seen that the trend becomes monotonic when particles are around 16 units large (Figure S7). See also Figure 5.

Analysis of the interactions in CHCl3

The systems under study present different hydroxylation patterns. In CHCl3, they may engage in extensive H-bonding during solvent evaporation and solid consolidation. While the conformational study supports that the -OH groups near the imine nitrogen atoms (i.e. position 2) are most likely involved in intramolecular interactions (in agreement with what has been observed in the crystal phase of related structures -Thakuria et al., 2012-), the ring decorated with 2 -OH (in 1 and 2) and 3 -OH moieties (in compound 3) in portion β of the scaffold is central to intermolecular H-bonding, giving place to multiple possible modes of aggregation relying on this directional interaction (Figure S8). Intermolecular H-bonding of 1, 2 and 3 is driven at the (sub)nanoscale regime by the accessibility of -OH groups in portion β and facilitated by the dielectric constant of the solvent, thus allowing multiple modes of aggregation during the earliest assembly stages, as shown in Figure S8 for molecule 2 in CHCl3. See also Figure 5.

While this kind of pattern is common to 2 and 3, due to the more accessible -OH group in position 4, that is not completely the case for compound 1. For this system, the outermost -OH is located at position 3, relatively hindered for H-bonding. Although the Turing pattern is not favoured in CHCl3, these results show that -OH at position 3 is still accessible for H-bonding, and this regulated hindrance (cf. -OH at position 3 in 1 vs. 4 in 2 and 3) seems to play an important role in the modulation of λ at the Turing patterns observed with each BB.

H-bond-driven aggregation can be safely assumed to appear since the earliest stages of particle formation, being a decisive factor during self-assembly, controlling the first hierarchy of supramolecular consolidation; this opens pathways for directional aggregation overcoming the isotropic scenario where diffusion and evaporation rate are probably the main factors controlling growth.

Analysis of the interactions in acetone

In acetone, intermolecular assembly through H-bonding can be expected to be only partially screened. Besides the difference in dielectric constant, the volatility of acetone, > 800 higher than that of DMSO, is likely an important factor during aggregation. Fast evaporation may lead mainly to kinetic ordering during aggregation (dynamic self-assembly). In contrast, the slow evaporation of DMSO may favour a thermodynamically driven assembly (static self-assembly). According to our findings, the Turing pattern is disfavoured in the latter conditions.

Once concentration reaches a certain value, screening of H-bonding may be overridden, opening access to the gamut of aggregation modes led by this interaction. As seen in Figure S9, aggregation of 1 occurs by partially hindered H-bonding domains; in the case, particle was modelled at the DFTB+ level with the 3ob Slater-Koster library, as detailed in the ‘in silico studies’ of this section. Regulation for H-bond-guided oligomerisation in molecule 1 appears as a crucial factor, in contrast with molecules 2 and 3 where functionalisation steric demand is poorer and the -OH group in position 4 in both BBs is more accessible to engage in H-bonded aggregates (Figure S10). In the left column of Figure S10, we find the most stable conformers for BBs 1 (top), 2 (centre) and 3 (bottom); In the right column of Figure S10, we find the respective modes of bimolecular aggregation by regulated access to intermolecular H-bonding for BBs 1 (top), 2 (centre) and 3 (bottom); in these cases, free molecules were modelled by MD in acetone and optimised at the DFT level and particles were modelled at the DFT level. Computational details can be found in the ‘in silico studies’ subheading of this section. See also Figure 5.

With these modelling experiments, we conclude that an intermediate dielectric constant and high volatility of the chosen solvent are the conditio sine qua non to achieve Turing-patterned self-assembled materials. From the molecular systems we explored, it is crucial that aggregation through H-bonding is modulated both sterically at the molecular unit and dynamically through judicious solvent selection.

Assembly dimensions in acetone

The largest particles we observed for 1–3, and whose assembly modelled computationally, lie around the 20–30 nm scale (Figures S11–S13). These dimensions are approximated by clusters of 210 molecular units, with, for instance, 418,251 uma for 1 (roughly >6 times the mass of the Human Serum Albumin protein); nonetheless their self-assembly and patterning via non-covalent interactions can be tuned through the adequate combination of physicochemical parameters.

In Figure S11, structures composed by one to 1024 molecules of 1 showing their solvent accessible surface are comparatively scaled; the number of assembled molecules is indicated below each agglomerate and their corresponding hierarchy levels are indicated above each group of structures. The structure containing one molecule conformer corresponds to the result obtained by MD in 1000 molecules of acetone, stripped from solvent molecules and further optimised at the DFT level. Clusters of up to 16 molecules were computed using DFTB+. All the supramolecular assemblies were optimised at the MM Dreiding level, using DFT-derived atomic charges, adjusting the number of optimisation steps to ensure minima were reached in all cases. Computational details are included in the ‘in silico studies’ subheading of this section. See also Figures 4 and 5.

In Figure S12, selected structures composed by one to 1024 molecules of 2 showing their solvent accessible surface are comparatively scaled; the number of assembled molecules is indicated below each agglomerate. Results were obtained with a similar procedure as in Figure S11. Computational details are included in the ‘in silico studies’ subheading. See also Figure 5.

In Figure S13, selected structures composed by one to 1024 molecules of 3 showing their solvent accessible surface are comparatively scaled; the number of assembled molecules is indicated below each agglomerate. Results were obtained with a similar procedure as in Figure S11. Computational details are included in the ‘in silico studies’ subheading. See also Figure 5.

Diffusivity approximations for the hierarchy level 1 (L1)

Diffusivity for L1 (upper limit, DL1 = 346 μm2s-1) was estimated using dL1 average from HR-TEM analysis data in Figure 4D, while its counterpart for L1 (lower limit, DL1 = 107 μm2s-1) was obtained using the estimated size d from molecular modelling size evolution trend for n = 256 in Figure 4E. Estimations were carried out using the Stokes-Einstein equation (Equation 1) at 298 K, making the following assumptions on the hydrodynamic diameter (dh) in acetone:

Supramolecular reaction mechanisms for 2-node RD network

Proposed supramolecular reaction mechanisms for the four interactions (Figure 6B) in the 2-node RD network between morphogens A and B (Figure S14), their kinetic constant predictions (κ) obtained from their estimated diffusivities DA = 17 μm2s-1 and DB = 1750 μm2s-1 and examples of experimental/molecular modelling/mathematical analysis data which might support the viability for the reaction steps proposed are listed ahead. In the mechanisms, level 1 is not accounted for as a diffusible participant in the overall RD network but is only limited to an intermediary role. The RD equations (Equations 2 and 3) derived from the automatic mathematical analysis are also listed below (Marcon et al., 2016).

In Figure S14, positive feedback (panel a top, green) from morphogen A (L2) to morphogen B (L0) might proceed in two reaction steps, i→i (bottom), where the overall reaction results into an increased local concentration of free molecules (L0). Negative feedback loop (panel a top, red) on morphogen B node could be carried out in one reaction a (bottom), which results into self-clearance of free BBs from the diffusible pool. Negative feedback (panel b top, red) from morphogen B to morphogen A potentially occurs in two reaction steps, i→i (bottom-right), which results into a decreased local concentration of L2 concerted with a local replenishment of BBs consumed in the process. Positive feedback loop (panel b top, green) on morphogen A node could proceed in three steps, i→a→a (bottom-left), where the overall reaction results into production of additional L2 particles at the diffusible pool. See also Figure 6.

κ. Positive feedback loop in morphogen A node (L2) through three successive reaction steps: i→a→a.

κ:0.187<κ<1 (×10−4 s−1).

i) Leaching of inhibited level 1 intermediaries (L1) from L2 without destruction of the latter.

a) Self-coalescence of intermediary L1 particles yields larger counterparts (activated level 1 intermediaries, L1).

a) Progressive assembly of intermediary L1 to generate L2 results into additional L2 particles at the diffusible pool.

Examples: Leaching of L1 can be observed during the bursting process in Figure 4C and their Gaussian trend in size is shown in Figure 4D. The homogeneous growth trend from L1 to L1 can be observed in molecular modelling, Figure 4E throughout the regression curve section between n = 8 and n = 256, whereas the isotropic assembling in the same n interval is evidenced in Figure S11. The assembly of level 1 structures to produce level 2 counterparts can be detected as a reassembly process in Figure 4C, while the regression curve trend disruption between n = 256 and n = 1024 in Figure 4E and anisotropic assembling in the same n interval in Figure S11 confirm the viability for such behaviour.

κ. Negative feedback from morphogen B (L0) to morphogen A in two reaction steps: i→i.

κ:–6.8<κ<−0.504 (×10−4 s−1).

i) BBs located at the interface of L2 assemblies could be accessible to sequestering by unassembled BBs (L0) to yield intermediary L1 nanostructures until remaining matter at L2 is consumed; in the process, destruction of L2 derives into a decreased local concentration of the same.

i) Previously formed L1 intermediaries that are unstable because of a low degree of oligomerisation could follow a splitting process to release unassembled molecules (L0) and thus, promote a local replenishment of BBs consumed at i.

Examples: The size reduction trend in L2 to values close to L1 observed for progressive dilution of BB in Figure 4A, as well as the bursting process detected in Figure 4C and the Gaussian trend in Figure 4D show that A is depleted by the local conditions through a process that is carried out via release of L1 particles to the surroundings. In this sense, in Figure 6D, the in-phase periodic drop in concentration of A and B (the troughs in the curves) where [A]local<<[B]local suggests that the higher availability of B triggers an interaction process between both morphogens which also drives the inhibition of A. In addition, the six scattered L1 structures in Figure 4C which appear markedly smaller than their reassembled counterparts, suggest that these L1 structures are potentially unstable if remained isolated from other L2 bodies (i.e. if the local concentration of L2 is relatively lower); in this sense, the less-marked concentration drop of B vs. A in Figure 6D, ([B]local>>[A]local) could be partially a result from a buffering effect from free molecules produced during the splitting of L1, which avoids a more noticeable depletion of B.

κ3: Positive feedback from morphogen A to morphogen B in two reaction steps: i→i.

κ:0.504<κ<6.8 (×10−4 s−1).

i) Leaching of L1 intermediaries from L2 assemblies without total destruction of the latter.

i) Successive splitting of unstable L1 into L0, should result in an overall increased local concentration of free molecules.

Examples: The bursting process of L2 in Figure 4C and the Gaussian trend in Figure 4D, together with the 6 scattered and potentially unstable L1 structures observed in Figure 4C suggest that L2 can release matter to the surroundings via smaller intermediary particles which potentially degrade into their elementary components. The overall process can be explained by the in-phase periodic rise in concentration between A and B in Figure 6D, where [A]local>>[B]local means that the higher the amount of A particles at a local level, the higher the local concentration of free molecules, therefore, these free molecules should be matter originated from L2 structures.

κ. Negative feedback loop in morphogen B node through one reaction: a.

κ:–1.98<κ<−0.504 (×10−4 s−1).

a) Progressive assembly of free molecules yields activated level 1 intermediaries (L1); the process results into self-clearance of free BBs from the diffusible pool.

Examples: In Figure 4E, the homogeneous trend observed throughout the regression curve section between n = 2 and n = 256 and the isotropic assembling in the same n interval in Figure S11 could explain the self-mediated clearance of free molecules.

Supramolecular reaction mechanisms for 3-node RD network

Proposed supramolecular reaction mechanisms for the six interactions (Figure 6C) in the 3-node RD network between morphogens A and B and hierarchy level 1 (Figure S15), their kinetic constant predictions (κ) obtained from the estimated diffusivities DA = 17 μm2s-1, DB = 1750 μm2s-1 and either L1 lower limit DL1 = 107 μm2s-1 or L1 upper limit DL1 = 346 μm2s-1, as well as examples of experimental/molecular modelling/mathematical analysis data which might support the viability for the reaction steps proposed are listed ahead. In the mechanisms, level 1 in either the activated L1 or the inhibited L1 form is accounted for as a diffusible and reactive component in the overall RD network, instead of an intermediary player as in the 2-node RD network equivalent. The RD equations (Equations 4, 5, and 6) derived from the automatic mathematical analysis are also listed below (Marcon et al., 2016).

In Figure S15, Positive feedback (a top row, green) from level 1 (either L1 or L1) to A (L2) might proceed through two possible routes (middle row), reaction a and/or reaction steps a→a, giving as a result an increased local concentration of A in both scenarios. Negative feedback (a top, red) from B (L0) to A could be carried out in two reaction steps i→i (middle row), which results into clearance of A from the diffusible pool concerted with a local replenishment of BBs consumed in the process. Positive feedback from A to level 1 (b top row, green) is potentially performed through three possible routes (middle row), reaction i, reactions i→a and/or reactions i→i→a; in all the cases, the local concentration of either L1 or L1 is increased at the diffusible pool. Negative feedback loop (b top row, red) on level 1 node could proceed in two possible reaction routes (bottom row), successive reaction steps i→i and/or reaction i; in both routes, self-promoted depletion of L1 and/or L1 might occur by instabilities arisen from local concentration conditions. Positive feedback (c top row, green) from level 1 to B is probably developed through three alternative routes (middle row), reaction steps i→i, only the reaction i and/or the reaction i; in the three scenarios, the local concentration of free BBs is increased by the presence of either L1 or L1. Negative feedback (c top row, red) from B to level 1 might be carried out through two possible routes, reaction steps i→i (bottom row) and/or only the reaction i; in both routes, free BBs are responsible for triggering the local depletion of any form of L1 (L1 or L1) from the diffusible pool. See also Figure 6.

κ. Negative feedback loop on hierarchy level 1 node (L1 and/or L1) via two possible routes occurring either individually or simultaneously: route departing from L1 in two-reaction steps i→i and route departing from L1 in one-reaction i.

For DL1 = 107 μm2 s−1, κ:–1.38<κ<−0.0173 (×10−4 s−1).

For DL1 = 346 μm2 s−1, κ:–51.8<κ<−0.5 (×10−4 s−1).

i) Local instabilities arisen from relatively low level 1 concentration conditions trigger splitting of L1 nanostructures to produce their smaller counterparts L1, these products participate in a further step i.

i) The low concentration of pre-existent and/or generated L1 bodies promotes their destabilisation which catalyses their splitting into free molecules to the surroundings. In the overall process, the concentration of both level 1 forms is decreased.

Examples: The six isolated level 1 nanostructures observed by TEM in Figure 4C appear smaller and less stable than their counterparts that remain relatively closer to other level 1 homologues. This suggests that the lower the concentration of level 1 nanostructures, the higher their probability to be degraded into the elementary components. The same phenomenon is supported by Figure 5C for self-assembly generations 2 and 3 (2 ≤ n ≤ 8), where their instability is predicted as a result from the energetic trade-off occurring between solvation and aggregation upon local BB dilution, while at this self-assembly generation range, H-bonding is the predominant motive force (Figures 5A, S9, and S10).

κ. Negative feedback from morphogen B (L0) to hierarchy level 1 (L1 and/or L1) via two possible routes occurring either individually or simultaneously: route starting from L1 in two-reaction steps i→i and route starting from L1 in one-reaction i.

For DL1 = 107 μm2 s−1, κ:–11.7<κ<−0.447 (×10−4 s−1).

For DL1 = 346 μm2 s−1, κ:–53.1<κ<−0.813 (×10−4 s−1).

i) BBs located at the interface of L1 nanostructures could be accessible to sequestering by unassembled molecules (L0) to yield oligomer counterparts L1 until remaining matter of L1 dimensions is consumed; the generated L1 particles participate in a further step i.

i) superficial molecules on pre-existent and/or generated L1 particles could be sequestered by free BBs to produce L1 counterparts at a lower degree of oligomerisation; these L1 products are susceptible to either progressive leaching of their components (releasing of L0) or to splitting into unassembled molecules due to their inherent instability, in both scenarios, the released L0 components may participate in further sequestering cycles with other L1 located at the surroundings. In the process, free molecules are responsible for triggering the destruction of both level 1 forms and thus, a decreased local concentration of the same, whereas in parallel, the L0 molecules consumed in each reaction are replenished to the medium.

Examples: The six isolated level 1 nanostructures observed by TEM in Figure 4C which are smaller than the average size measured in Figure 4D indicate that they are potentially unstable because of their reduced local concentration. However, the coexistence of free molecules in solution (e.g. free molecules detected by DOSY at C in Figure S2) together with the in-phase drop in concentration between level 1 and B observed in Figure 6E (for L1) and Figure S16A (for L1), where [B]local>>[L1]local or [B]local>[L1]local, suggest that the relatively higher concentration of morphogen B should catalyse local depletion of level 1 through physical interactions between both components. According to molecular modelling of L1 particles with a low degree of oligomerisation (Figures 5A, S10, and S11 and our analysis of the interactions in acetone), it is expected that throughout the sequestering process from their larger counterparts, H-bonding played a similar leading role as in the stabilisation of L1 during its formation via self-assembly. Moreover, according to the solvation vs. aggregation energy profile in Figure 5C, the rise in the aggregation energy component for the L1→L1 range (towards less exothermic values) is less marked than the drop in the solvation energy component at the same interval and direction (towards more exothermic values). This suggests that at these regimes, highly solvated/more concentrated free molecules could participate with a higher probability and more actively in interfacial sequestering to yield smaller oligomers which gain on solvation energy component at the expense of parent L1 particles.

κ. Positive feedback from hierarchy level 1 (L1 and/or L1) to morphogen B via three alternative routes occurring either individually or simultaneously: two-reaction route i→i, one-reaction route i and one-reaction route i.

For DL1 = 107 μm2 s−1, κ:0.00173<κ<0.00904 (×10−4 s−1).

For DL1 = 346 μm2 s−1, κ:0.769<κ<12.8 (×10−4 s−1).

i) if the local concentration of BBs in the form of L1 surpasses the local concentration of L0, the level 1 nanostructures exposed to these conditions could spontaneously carry out a leaching process where free molecules are released to the surroundings, and smaller L1 counterparts are obtained as a by-product; the latter may or may not participate in a further step i.

i) if the local concentration of pre-existent and/or generated L1 exceeds the concentration of L0, L1 particles could go on with the leaching process to produce L1 counterparts at a progressively lower degree of oligomerisation. In the process, additional free molecules are constantly released to the medium (increase in L0 concentration) by the presence of level 1 structures.

Examples: The in-phase periodic rise in concentration of level 1 and morphogen B at a local level in Figure 6E (L1) and Figure S16A (L1) where [L1]local>>[B]local or [L1]local>[B]local, denotes that the increased concentration of morphogen B is determined by a higher availability of level 1 structures which give birth to the first ones; in such process, a progressive leaching mechanism is the most probable route. This scenario is supported by Figure 5C, where instability in self-assembly arises in a natural manner as long as the degree of oligomerisation decreases for both L1 forms, due to an energetic trade-off between solvation and aggregation which predominantly yields L0.

κ. Positive feedback from morphogen A (L2) to hierarchy level 1 (L1 and/or L1) via three alternative routes occurring either individually or simultaneously: one-reaction route i, two-reaction route i→a, and three-reaction route i→i→a.

For DL1 = 107 μm2 s−1, κ:0.712<κ<2.24 (×10−4 s−1).

For DL1 = 346 μm2 s−1, κ:0.0479<κ<1.23 (×10−4 s−1).

i) generation of L1 nanostructures originated from kinetically unstable L2 particles could spontaneously be carried out through a leaching process where level 1 particles are released to the surroundings, as long as precursor L2 particles gain enough local stability to maintain their hierarchy level, so as to avoid a total degradation of the same; depending on the local conditions, leached L1 may or may not participate in either i or a as a further step, however, in this step, a higher concentration of level 1 (inhibited form) is triggered by the only presence of level 2.

a) if the local conditions are favourable, the previously leaked L1 structures may proceed with a spontaneous self-coalescence to yield level 1 counterparts with a higher oligomerisation degree (L1). As an overall result, concentration of level 1 (activated form) is also risen by L2.

i) if L1 particles expelled through i are exposed to locally unfavourable conditions that promote their kinetic instability, they may spontaneously split into their components; this process promotes a transient increase in free molecules which participate in a further step a.

a) the transient rise in concentration of free molecules produced during i under unstable conditions favours their rapid (and dynamic) assembly into L1 particles. The global result of this route is also an increased concentration of level 1 (activated form) started by L2.

Examples: The bursting process observed in Figure 4C indicates that L1 particles may be produced from L2 precursors without destruction of the latter. Either the assembly of free molecules or the coalescence of L1 structures to yield higher degrees of oligomerisation (e.g. L1) is viable according to molecular modelling results; such process can be carried out starting from an oligomerisation degree as low as n = 2 (Figures 4E, 5A, S9, and S11) to rapidly produce more-stable particles with n ≥ 8 (Figure 5C). This process is driven at the first stages by H-bonding between portion β in 1, but, as oligomerisation increases, vdW and π-π interactions through portion α are increasingly playing as the leading motive force. The in-phase periodic rise in concentration of level 1 and morphogen A at a local level in Figure 6E (L1) and Figure S16A (L1) where [A]local>[L1]local, indicates an activation of L1 by the higher concentration of morphogen A.

κ. Positive feedback from hierarchy level 1 (L1 and/or L1) to morphogen A through two possible routes taking place either individually or simultaneously: two-reaction route a→a and one-reaction route a.

For DL1 = 107 μm2 s−1, κ:0.00139<κ<0.00387 (×10−4 s−1).

For DL1 = 346 μm2 s−1, κ:0.119<κ<0.813 (×10−4 s−1).

a) depending on the local conditions, coalescence of L1 particles could produce level 1 versions in progressively higher oligomerisation degrees at L1. These local conditions might be favoured at high L1 concentration values. The L1 products participate in a further step a.

a) L1 nanostructures which are either pre-existent or produced in a previous step a, could generate L2 particles if their local concentration is high enough so as to facilitate their successive assembly. In the process, L2 concentration is increased by level 1 hierarchy whether in its inhibited or its activated version.

Examples: The reassembly process observed in Figure 4C indicates that free level 1 nanostructures are responsible for the formation of L2 bodies as long as their local concentration is high enough so as to trigger the process (e.g. their spatial distribution is constrained to a reduced region). Formation of L2 from level 1 nanostructures and the coalescence of L1 into L1 are also demonstrated by molecular modelling (Figures 4E, 5B, and S11), where π-π and vdW interactions through portion α in 1 play a leading role (agonistic signalling). Either activated or inhibited level 1 structures carry out an activation task on morphogen A, according to numerical simulation results in Figure 6E (for L1) and Figure S16A (for L1); in these results, the in-phase periodic oscillations in both concentrations show that a locally increased amount of level 1 sharply triggers a higher concentration of level 2. The overall process is synergistically connected to the opposite direction feedback, where A activates level 1 (κ).

κ. Negative feedback from morphogen B to morphogen A through a two-reaction route: i→i.

For DL1 = 107 μm2 s−1, κ:–0.018305<κ (×10−4 s−1).

For DL1 = 346 μm2 s−1, κ:–0.375<κ (×10−4 s−1).

i) BBs located at the interface of L2 particles could be accessible to sequestering by unassembled molecules (L0) to yield oligomer counterparts L1. This process continues until remaining matter of L2 dimensions is consumed and only L1 are present at a local level; in the process, destruction of L2 derives into a decreased local concentration of the same; the generated L1 particles participate in a further step i.

i) the local concentration conditions and low degree of oligomerisation of the previously formed L1 structures facilitate their destabilisation and further splitting into free molecules which are released to the surroundings. In the process, L0 molecules consumed during the first step are replenished to the diffusible pool, whereas local concentration of L2 is depleted as a result of its interaction with free molecules.

Examples: As in the mechanism for κ at the 2-node RD network scenario, the size reduction trend in L2 to values close to L1 observed for progressive dilutions in Figure 4A, in addition to the bursting process observed in Figure 4C and the Gaussian trend in Figure 4D indicate that A is depleted via release of L1 particles to the environment. In this regard, since L1 structures possess a low degree of oligomerisation, their transient stabilisation is mainly led by H-bonding through portion β in 1, according to molecular modelling (Figures 5A and S11 and the analysis of the interactions in acetone), therefore, it is expected that the sequestering process might be also driven by this non-covalent interaction (antagonistic signalling). However, as long as the produced L1 particles remain at low local concentrations, the trade-off occurring in solvation and aggregation energetics (Figure 5C) compromises their stability, thus displacing equilibria to their splitting product L0. Similarly, the analysis on the solvation vs. aggregation energy curves in Figure 5C also shows that the increase in aggregation energy for the L2→L1 range (to less exothermic values) is not as steep as the drop in solvation energy at the same range and direction (to more exothermic values). This suggests that upon local dilution regimes, highly solvated/more concentrated free molecules could participate with a higher probability and more actively in interfacial sequestering to yield smaller oligomers which gain on solvation energy component at the expense of; parent L2 particles. As shown in Figures 6E and S16A, the in-phase periodic drop in concentration of A and B where [A]local<<[B]local suggests that the higher availability of B triggers an interaction process between both morphogens which also drives the inhibition of A. In addition, the smaller size of the six L1 structures in Figure 4C (which appear isolated from the rest of particles), suggests that these bodies are potentially unstable if remain distant from other L1 or L2 bodies (i.e. if the local concentration of L1 or L2 is relatively lower); this situation can be explained by the in-phase drop in concentration of all the diffusible components A, B and L1 in Figures 6E and S16A, where the ratio [B]local>[L1]local>[A]local could indicate that free molecules produced during the splitting of L1, avoid a scenario with a more marked depletion of B.

Solvation energies on supramolecular clusters

We were able to estimate semi-quantitative solvation energies with respect to the SMD method by leveraging on the Dreiding force field for geometry optimisations and energy calculations. ESP-mapped atomic charges obtained by DFT with the PBE(DTS)/TNP with COSMO solvation proved to yield the best accuracy and were used throughout all our studies. We found that the molecular silverware procedure (MSP) provides excellent results when combined with the selected forcefield and coordination numbers are adequately converged.

The method for calculation of solvation energies is performed as follows:

Generated optimised geometries and atomic charges for both the solute and the solvent, labelled as ‘screen’ and ‘base’ in the procedure, respectively.

The number of iterations (N) in the MSP step depends on the relative sizes of the Base and the Screen and must be determined for each generation. As an example, we show the sweep from N = 30 to N = 350. We emphasise that extensive tests proved that three runs at N = 30, 50 and 150 represent an adequate compromise, saving computation time while providing coordination numbers within 0.2 units from more exhaustive sampling tests. As seen in Table S1 and Figure S21, the coordination number (C) at infinite N is obtained by simple extrapolation. In that example, the extrapolated C = 22.5.

Three mixing runs must be performed with rigorous sampling. We propose N = 30, with 106 energy samples and an energy bin width of 0.02 kcal/mol with the selected FF. These parameters were adequate for the entire set of particles studied but must be critically selected based on the deviation of individual results from the average. This step yields the average energy value for Base + Screen interaction, EBS.

Solvation energy is obtained as ΔEsolv = EBS × CN.

For generations 8, 9 and 10 (256, 512 and 1024 molecules), we leveraged on the linear relationship between aggregate size and the weighted solvation energy as observed in Figure 5C. We thus performed a linear extrapolation as shown in Figure S22 that afforded information about the earlier generations. Weighted solvation energy is inversely proportional to the number of molecules per aggregate at earlier generations. This behaviour is useful to assess the system at later generations and larger aggregate sizes.

We estimated the weighted solvation energies for generations 8, 9 and 10, represented in Figure 5C as hollow triangles, by applying this simple linear model. All solvation energies obtained are condensed in Table S2, expressed in kcal/mol.

Aggregation energies on supramolecular clusters

Aggregation energies were obtained through the molecular silverware methodology. Due to the additivity of this protocol, total aggregation energies were calculated additively for all generations (number of molecules, n > 1) by application of the general Equation 7:

In Equation 7, n is the target generation, m are lower generations, N is the number of molecules for a given generation and E is the molecular silverware energy term. Naturally, there is no aggregation energy for generation 0, whereas for generation 1 it is simply the Base + Screen interaction energy . Thus, this formula is useful only for higher generations. The corresponding energy values for all generations are in Table S3.

From these values, the total and weighted aggregation energies can be derived straightforwardly and are shown in Table S4 and Figure 5. All energies in both tables are expressed in kcal/mol.

Quantification and statistical analysis

Statistical analysis of micrographic data from OM, FFT, CLSM and HRTEM, was performed by direct inspection and manual collection of points and/or length values and/or graphical thresholds and/or maximum intensity trajectories using the ImageJ tools Measure, Analyze Particles, Histogram, Distribution, and Plot Profile, as provided within the 1.50b version; graphical thresholds in the annular spots from FFT represented the rim of the ring’s edges located at both, the internal and the external circumferences, where the intensity values dropped to a minimum equivalent to the background baseline value; maximum intensity trajectories in the annular spots from FFT typically represented the midpoint circumference located between the two graphical thresholds mentioned earlier. Statistical analysis of surface profile data from AFM was carried out by manual collection of distances between crests from data achieved with XEI package at the sampled trajectories. Average values reported in Figures 3G, 3I, and 4D, are arithmetic means of a set of values in a sampling number (N); central values indicated for FFT in Figures 1B, 3D, 6F, 6G, S1F, S1G, and S16B, are the maximum intensity trajectories in the annular spots as explained before; tolerance values in the figures are either range limits dictated by the graphical thresholds (Figures 1B, 3D, 6F, 6G, S1F, S1G, and S16B), standard deviation - s.d.- of a set of values in a sampling number -N- (Figures 3I and 4D), or confidence intervals (CI) at the indicated degrees of freedom (ν) for N = ν+1, under a t-distribution scheme (Figure 3G), as pointed out in each panel. CI values are reported for p < 0.02, as described accordingly. Gaussian curve approximations in Figures 3I and 4D represent the normal distribution fit for the shown histograms at the found arithmetic mean values, s.d. and N, as indicated in the respective figure panels. Size of N was chosen by visual inspection and represented the total number of collectable elements in the sample (Figures 3G, 3I, and 4D) or the number of collectable points that rendered a constant arithmetic mean value and simultaneously met the respective graphical criteria mentioned earlier (Figure 3D).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Chemicals, peptides, and recombinant proteins
2-hidroxybenzophenone	Sigma-Aldrich	CAS: 117-99-7
2-hydroxybenzaldehyde	Sigma-Aldrich	CAS: 90-02-8
Methanol	Sigma-Aldrich	CAS: 67-56-1
o-phenylenediamine	Sigma-Aldrich	CAS: 95-54-5
2,3-dihydroxybenzaldehyde	Sigma-Aldrich	CAS: 24677-78-9
Acetone-d6	Sigma-Aldrich	CAS: 666-52-4
CDCl3	Sigma-Aldrich	CAS: 865-49-6
2,4-dihydroxybenzaldehyde	Sigma-Aldrich	CAS: 95-01-2
2,3,4-trihydroxybenzaldehyde	Sigma-Aldrich	CAS: 2144-08-3
2,5-dihydroxybenzaldehyde	Sigma-Aldrich	CAS: 1194-98-5
Acetone	Sigma-Aldrich	CAS: 67-64-1
Chloroform	Sigma-Aldrich	CAS: 67-66-3
Dimethyl sulfoxide	Sigma-Aldrich	CAS: 67-68-5
3-((E)-((2-(((E)-(2-hydroxyphenylphenyl)methylene)amino)phenyl)imino)methyl)benzene-1,2-diol	This paper	N/A
4-((E)-((2-(((E)-(2-hydroxyphenyl) (phenyl)methylene)amino)phenyl)imino)methyl)benzene-1,3-diol	This paper	N/A
4-((E)-((2-(((E)-(2-hydroxyphenyl) (phenyl)methylene)amino)phenyl)imino)methyl)benzene-1,2,3-triol	This paper	N/A
2-((E)-((2-(((E)-(2-hydroxyphenyl) (phenyl)methylene)amino)phenyl)imino)methyl)benzene-1,4-diol	This paper	N/A
Software and algorithms
ImageJ 1.50b	Wayne Rasband	https://imagej.nih.gov/ij/index.html
FV10-ASW Ver. 01.07c	Olympus Corporation	https://www.olympus-lifescience.com/es/support/downloads/
XEI Package	Park Systems Corporation	https://www.parksystems.com/manuals-software
Gatan Digital Micrograph	Gatan, Inc.	https://www.gatan.com/products/tem-analysis/digitalmicrograph-software
Materials Studio 8 suite	Biovia, Dassault Systèmes	https://www.3ds.com/products-services/biovia/products/molecular-modeling-simulation/biovia-materials-studio/
ORCA 4.0 program	Neese (2012)	https://wires.onlinelibrary.wiley.com/doi/10.1002/wcms.81
Molecular Silverware	Blanco (1991)	https://onlinelibrary.wiley.com/doi/10.1002/jcc.540120214
RDNets	Marcon et al. (2016)	https://elifesciences.org/articles/14022
Wolfram CDF Player 11.3	Wolfram Research, Inc.	https://www.wolfram.com/cdf-player/index.es.html
Other
Perkin Elmer 2400 for CHNS analyser	Perkin Elmer, Inc.	https://www.perkinelmer.com/es/product/2400-chns-o-series-ii-system-100v-n2410650
Varian VNMRS 400 MHz apparatus	Varian Instruments, Inc.	https://www.varian.com/
Agilent 6530 QTOF spectrometer	Agilent Technologies, Inc.	https://www.agilent.com/en/product/liquid-chromatography-mass-spectrometry-lc-ms/lc-ms-instruments/quadrupole-time-of-flight-lc-ms/6530-q-tof-lc-ms
Bruker Ascend 500 MHz apparatus	Bruker Corporation	https://www.bruker.com/en/products-and-solutions/mr/nmr/ascend-nmr-magnets.html
Motic model BA210	Motic	https://www.motic.com/As_LifeSciences_UM_BA210/
Olympus model FV1000 microscope	Olympus Corporation	https://www.olympus-lifescience.com/es/technology/museum/micro/2004/
Park Systems model NX10 microscope	Park Systems Corporation	https://www.parksystems.com/products/small-sample-afm/park-nx10/specifications
Zetasizer Nano ZSP	Malvern Instruments Ltd.	https://www.malvernpanalytical.com/es/support/product-support/zetasizer-range/zetasizer-nano-range/zetasizer-nano-zsp
JEOL model JEM-2010 microscope	JEOL, Ltd.	https://www.jeol.co.jp/en/products/