A stochastic model study on the self-assembly process of a Pd2L4 cage consisting of rigid ditopic ligands

Satoshi Takahashi *a, Yuya Sasaki a, Shuichi Hiraoka *a and Hirofumi Sato *bc
aDepartment of Basic Science, Graduate School of Arts and Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan. E-mail: takahasi@mns2.c.u-tokyo.ac.jp; chiraoka@mail.ecc.u-tokyo.ac.jp
bDepartment of Molecular Engineering, Kyoto University, Kyoto 615-8510, Japan. E-mail: hirofumi@moleng.kyoto-u.ac.jp
cElements Strategy Initiative for Catalysts and Batteries, Kyoto University, Kyoto 615-8510, Japan

Received 28th September 2018 , Accepted 28th November 2018

First published on 28th November 2018

The coordination self-assembly process of a Pd2L4 cage including rigid ditopic ligands, L, was studied numerically. A recently developed experimental approach (QASAP: quantitative analysis of self-assembly process) revealed that the rate-determining steps in the self-assembly of the Pd2L4 cage are intramolecular ligand exchanges at the late stages of the self-assembly. In this study, the self-assembly process before the rate-determining steps, which could not be investigated by experiment, was analyzed based on a minimal reaction network model. Only eight variable parameters of rate constants for ligand exchange reactions are sufficient enough to reproduce the time evolution of substrates and the products during the self-assembly of the cage. With these parameters, the major self-assembly pathway was determined. It was also found that a non-negligible amount of an incomplete cage (IC), Pd2L3X2 (X indicates the leaving ligand), which was not suggested by QASAP, should be transiently produced. Numerical tests also suggest that the small rate constant value of the intramolecular ligand exchanges due to a restricted geometry causes the final stage to seemingly become the rate-determining step.

1 Introduction

Self-assembly of molecular building blocks is one of the most important reaction processes to construct 3-dimensional nanostructures. One of the most reliable ways to reveal molecular self-assembly processes is to obtain information about all the intermediates produced during the self-assembly. However, it is practically impossible to monitor the time evolution of all the species transiently produced in the self-assembly. In such situations, Hiraoka and coworkers have recently developed a method that enables one to investigate coordination self-assembly processes by monitoring the time evolution of the average composition of all the reaction intermediates (QASAP, quantitative analysis of self-assembly process).1–3 This method tracks self-assembly by the quantification of all the observable chemical species (mainly reactants and products) instead of the detection of the intermediates. The information obtained from reactants, products, and long-lived intermediates allows one to monitor a time-dependent concentration of non-observable species by the principles of mass balance. Coordination self-assemblies whose reaction processes have been investigated by QASAP so far include Pd(II)-linked cages,4–8 octahedron-shaped capsules,9 a tetrahedron,10 a sphere,11 a double walled triangle12 and a square,13 and Pt(II)-linked hexagonal macrocycles.14,15 Among them, the Pd2L4 cages composed of a total of six components are the simplest self-assembly systems. Despite the simplicity of the Pd2L4 cages, their formation processes are not so simple as to be described by a single reaction pathway.4 Indeed, QASAP revealed the self-assembly processes of Pd2L4 cages, but tracking all the intermediates is impossible not only by QASAP but also by any other experimental techniques. Thus, theoretical approaches are inevitable to fully understand molecular self-assembly processes.

In this study, based on a stochastic approach using the chemical master equation, the self-assembly process was revealed numerically through constructing a reaction model for the Pd2L4 cage formation. By using the stochastic simulation, it becomes possible to track the reactions whose time scales are much shorter than those monitored by 1H NMR spectroscopy. Details of the dynamics of reaction intermediates obtained by QASAP are further analyzed by investigating the time evolution of individual molecular species. Here, we applied a stochastic approach to the Pd214 system (Fig. 1), in which a rigid ditopic ligand 1 is composed of two 3-pyridyl groups connected to a benzene ring through ethynylene spacers. In this cage formation reaction the rate-determining step was experimentally identified clearly as the final stages of the self-assembly process. Let us consider the self-assembly of the Pd214 cage expressed by the following equation,

2·[PdX4]2+ + 4·1 ⇄ [Pd214]4+ + 8·X,(1)
where X is the leaving ligand that originally coordinates to the Pd(II) center of the metal source. The self-assembly of the Pd214 cage takes place through intra- and intermolecular multiple ligand exchanges of X with a ditopic ligand 1. Here we study the self-assembly process of the Pd214 cage from image file: c8cp06102e-t1.tif (Py* indicates 3-chloropyridine, which was used as a leaving ligand X in eqn (1)) and rigid ditopic ligand 1 (Fig. 1(a)).

image file: c8cp06102e-f1.tif
Fig. 1 Self-assembly process of the Pd2L4 cage. (a) Reaction formula and component molecules. In this numerical analysis, rigid ditopic ligand 1, in which two 3-pyridyl groups are connected to a benzene ring through ethynylene spacers, were chosen as L, and the leaving ligand X is 3-chloropyridine. (b) A schematic representation of reaction pathways revealed with our stochastic simulations. Thick black arrows indicate the main reaction route from the reactants to the products, while thin gray ones correspond to the minor pathways present in a non-negligible amount.

All the species in the self-assembly of the Pd214 cage should be described as Pda1bXc (a, b, and c are the numbers of components and thus have integer values), which is also described as (a,b,c) in this paper. Although it is impossible to follow the time evolution of all the intermediate species, 1H NMR spectroscopy makes it possible to quantify all the reactants and products, which in turn enables one to monitor the average composition of the intermediates, Pda1bXc. To discuss the self-assembly process the following two parameters n and k are defined by the following equations:

image file: c8cp06102e-t2.tif(2)
image file: c8cp06102e-t3.tif(3)
N is the coordination number of the metal ion and N = 4 in this study. While the parameter n represents the average number of metal ions bound to the ditopic ligand in Pda1bXc (in the present case, the maximum value of n is 2, the number of coordination sites of the ditopic ligand L), k is the average ratio of the metal ions to the ditopic ligands in Pda1bXc. As to the average composition of all the intermediates, Pda1bXc, which can be determined by QASAP, the 〈n〉 and 〈k〉 values are expressed in the same way,
image file: c8cp06102e-t4.tif(4)
image file: c8cp06102e-t5.tif(5)
The time evolution of the (〈n〉, 〈k〉) values on the (n–k) map provides us with useful information on the self-assembly processes, for which we attempt to obtain further details by employing the numerical approach.

The remaining part of the present article is organized as follows. Section 2 presents the procedure to perform the numerical calculations based on the stochastic approach. Numerical results are shown in Section 3, with comparisons to the original experimental data. Our numerical approach revealed early stages of the self-assembly process of the Pd214 cage that could not be investigated by QASAP. Finally, Section 4 gives the discussion and concludes this paper.

2 Preparation for numerical calculations

In our minimal network model, eight variable reaction rate constants (k1,…, k−4) were defined for three types of typical ligand exchanges and the final step was treated separately (see the minimal reaction network in Fig. 2).
image file: c8cp06102e-f2.tif
Fig. 2 Minimal reaction network in our stochastic simulations. Yellow, blue, and green pictorial figures represent the Pd(II) center, the ditopic ligand L, and the leaving ligand X, respectively. Red, blue, green, and yellow lines correspond to the intermolecular ligand exchange, bond formation of two Pd(II) centers, the intramolecular ligand exchange, and the final step leading to the cage, respectively. Each species is given a number for the sake of references in the main text.

• For the intermolecular ligand exchange reactions, [PdaLbXc]2a+ + L → [PdaLb+1Xc−1]2a+ + X, a = 1 or 2: k1 [min−1 M−1] and k−1 [min−1 M−1] for forward and backward reactions, respectively.

• For the formation of dinuclear species from two molecules of mononuclear species, whose Pd(II) ions are bridged by ditopic ligand L, [PdLb1Xc1]2+ + [PdLb2Xc2]2+ → [Pd2Lb1+b2Xc1+c2−1]4+ + X: k2 [min−1 M−1] and k−2 [min−1 M−1] for forward and backward reactions, respectively.

• For the intramolecular ligand exchange reactions in dinuclear species, [Pd2LbXc]4+ → [Pd2LbXc−1]4+ + X: k3 [min−1] and k−3 [min−1 M−1] for forward and backward reactions, respectively.

• Rate constants for the reaction [Pd214X]4+ → [Pd214]4+ + X are handled separately due to the experimental fact that this step was found to be rate-determining: k4 [min−1] and k−4 [min−1 M−1] for forward and backward reactions, respectively.

Note here that each rate constant is defined per reaction site, based on the modeling procedure reported previously.16 So the actual reaction rate for each reaction is estimated as the above constant multiplied by the total number of available combinations. For example, for a ligand exchange reaction between [PdX4]2+ (#1 in Fig. 2) and 1 to produce [Pd1X3]2+ (#2) and X, the rate constant is given as k1 times 4 (the number of Pd(II)–X bonds in [PdX4]2+) times 2 (the number of free pyridyl nitrogen atoms in 1), i.e.,

image file: c8cp06102e-t6.tif(6)
We adopted this setting to explicitly distinguish the structural difference among the species with the same composition. According to the experimental result that the self-assembly of the Pd214 cage mainly takes place through intermediates that do not contain more component than the cage (a ≤ 2, b ≤ 4),4 we have omitted the possibility of containing intermediates with three and more Pd(II) centers in this model. It is true that this network model simplifies the real system, but as we discuss later, even this simplified procedure gives valuable information about the species produced in a very early stage of the self-assembly, which should be a good starting point for the discussion of the essence of the molecular self-assembly process.

We constructed one of the simplest reaction networks as follows: Starting from the Pd214 cage (2,4,0), the reaction path is traced back to the reactants, i.e., [PdX4]2+ and 1. Those intermediates inaccessible “directly” (that is, via a single step) in the course of this tracing-back were not considered as the intermediate species responsible for the cage formation. For example, while (2,4,1) is included in the reaction network, (2,5,1) is not considered as a participant due to its inaccessibility in-one-step from (2,4,0), although both of them are accessible from (2,4,2). This procedure agrees with the consideration of the original article in that only four (#21, #22, #23 and #24 in Fig. 2) out of six (2,4,2) conformational isomers were considered as potential candidates for the precursor of the final cage. With a free rotation around the Pd–1 bonds assumed, and with the enantiomers treated as the same species, the number of species adopted for the overall reaction is 29, including the reactants and products. It was found that this network contains the total number of 136 reactions. A schematic representation of the reaction network is given in Fig. 2.

In order to track the time evolution numerically, there are some numerical approaches available. As performed in the original article, it is possible to directly follow the time evolution of the concentration of each species with a set of ordinary differential equations.4 However, to avoid the numerical instability caused by significantly different time scales expected to appear in different kinds of reactions, we have adopted the chemical master equation approach, the so-called Gillespie algorithm.17–20 In this algorithm, for all the possible M chemical reactions including molecular species Sai, Sbi, Sci, …,

Sai + Sbi + ⋯ → Sci + ⋯ (i = 1,…, M),(7)
the total reaction rate Rtot is calculated as
Rtot = r1 + r2 + ⋯ + rM,(8)
ri = ki[Sai][Sbi]….(9)
Starting from the initial time t = 0, at each instant t, which one of the reactions occurs is determined with the uniform random number s1 ∈ (0,1) as
image file: c8cp06102e-t7.tif(10)
Another uniform random number s2 ∈ (0,1) is independently given to fix the time incremental dt as
dt = ln(1/s2)/Rtot.(11)
Time is updated as t = t + dt, together with the update of the numbers of corresponding molecular species, i.e., 〈Sai〉 → 〈Sai〉 − 1, 〈Sbi〉 → 〈Sbi〉 − 1, 〈Sci〉 → 〈Sci〉 + 1, …. The reason why this algorithm traces the chemical reactions and actually works well is given in the literature in detail.17 This stochastic algorithm was previously applied to the self-assembly process of Pd(II)-linked octahedron-shaped capsules and successfully reproduced the experimental data by considering all the intermediate species of different compositions.21 In this study, we considered possible structural isomers with the same composition that can finally lead to the Pd2L4 cage.

At t = 0, initial values are given only to two components, that is, 〈[PdX4]2+0 = 1275 and 〈10 = 2550. These values are calculated from the averaged experimental initial concentrations for both components, [Pd2+]0 = 0.853 mM and [1]0 = 1.71 mM.4 Avogadro's constant is set to be NA = 6.0 × 1023 mol−1. Accordingly the volume of the simulation box was taken to be 2.5 × 10−18 L. The order of the initial numbers of molecules is determined based on the study of the ligand exchange model, in which it was numerically confirmed that ∼103 particles are enough to fix the rate constant values.16

3 Results

3.1 Reproduction of the experimental data with the minimal reaction network

First, we show the numerical results obtained from the minimal reaction network model. For this reaction network, by fitting the numerical data of (i) the existence ratios of all the reactants and products and (ii) 〈n〉 and 〈k〉 values along the time evolution to the experimental counterparts, rate constant searches in the parameter space were performed in order to reproduce the experimental results reasonably well. The adequacy of the fitting was evaluated from the residual sum of squares (RSS) with the experimental data of three runs (see the ESI of ref. 4). That is, for all the time steps ti at which the experimental existence ratio data RSex and parameters 〈nex and 〈kex are available, RSSs are calculated with the numerically obtained values RSnu as (note that S = [PdX4]2+, 1, [Pd2L4]4+, or X)
image file: c8cp06102e-t8.tif(12)
image file: c8cp06102e-t9.tif(13)

Fig. 3(a and b) shows examples of the time evolution of the existence ratios of [Pd214]4+ and X, and that of 〈n〉, respectively, both obtained numerically. Here, the initial numbers of molecules were set to be 10 times larger in order to make the fluctuations as small as possible, that is, 〈[PdX4]2+0 = 12[thin space (1/6-em)]750 and 〈10 = 25[thin space (1/6-em)]500. In Fig. 3(a), it can easily be seen that the behaviors of the existence ratio for both species are well reproduced by using only eight rate constant values. Note that the existence ratios of reactants were also confirmed to give good agreement with the experimental results, although they are not shown in the same figure due to their too simple behaviors (i.e., too rapid and complete depletion within the initial 5 minutes). For this result to well agree with the experimental counterpart (note that only the average data of three experimental runs with error bars is plotted for comparison), rate constants were set to be (k1, k−1, k2, k−2, k3, k−3, k4, k−4) = (103.6, 10−0.1, 104.0, 10−0.4, 10−1.0, 10−2.0, 10−0.9, 10−0.3). This set of rate constant values indicates the rapid ligand replacement of the ditopic ligand L with the leaving ligand X. Using the same set of rate constants, Fig. 3(b) shows good agreement of the time evolution, especially the initial rapid rising, of 〈n〉. The corresponding nk plots are shown in Fig. 3(c). It should be noted that the 〈k〉 value is also reproduced well with this rate constant set.

image file: c8cp06102e-f3.tif
Fig. 3 Fittings of numerical results to the experimental counterparts. Numerically determined rate constants are (k1, k−1, k2, k−2, k3, k−3, k4, k−4) = (103.6, 10−0.1, 104.0, 10−0.4, 10−1.0, 10−2.0, 10−0.9, 10−0.3). (a) Time evolution of the existence ratios for the product species. Blue and red curves correspond to [Pd214]4+ and X, respectively. Discrete plots (green, ×) are the experimental counterparts. (b) Time evolution of 〈n〉 obtained experimentally (green, ×) and numerically (red, +). Note that the sharp initial rise during the first ∼10 minutes is well reproduced with only eight parameters. (c) Numerically obtained nk plot (red, +). The experimental result is also plotted with green ×-marks for comparison (average of three runs in ref. 4 with error bars). Blue dashed circles indicate the locations corresponding to the (n,k) values for (2,4,2) and (2,4,1). It should be emphasized that even in the presence of the (2,3,2) incomplete cage (IC) the k value does not deviate from k = 0.5 (see also the main text).

3.2 Time evolution of individual molecular species

Next, let us inspect the details of the self-assembly process, by tracing time evolutions of individual components. These are shown in Fig. 4 with the logarithm of time as the abscissa axis. From this figure, we classify the self-assembly process of [Pd214]2+ cages into the following three stages:
image file: c8cp06102e-f4.tif
Fig. 4 Time evolution of the numbers of individual species. (a) Reactants and mononuclear species. Note that the numbers of (1,0,4) and 1 come down from 〈[PdX4]2+0 = 12[thin space (1/6-em)]750 and 〈L〉0 = 25[thin space (1/6-em)]500, respectively. (b) Dinuclear species (2,1,6) and (2,2,5) with only a single bridge. (c) Due to the large rate constant value k1 of the intermolecular ligand exchange, intermediate species with a single bridge (2,3,4) and (2,4,3) increase faster. (d) After almost complete consumption of the reactants, the formation of intermediates with the second and the third bridge follows. (e) The final Pd214 cage.

• In the early stage, a mononuclear species (1,1,3) with a single ditopic ligand 1 is formed (#2), followed by the rapid formation of binuclear species (2,1,6) and (2,2,5) (#7, #8, and #9) with only a single bridge (<10−1 min). Small numbers of the other mononuclear species (1, b, 4 − b), 1 < b ≤ 4 (#3, #4, and #5) are derived from the fact that k1 is smaller than k2 (Fig. 4(a and b)).

• Dinuclear species react with the free ditopic ligand 1, leading to almost complete depletion of them (10−1–100 min). In this process, the number of bridges is kept to one. Subsequent intramolecular ligand exchanges slowly give precursors of the cage (#21, #22, #23, and #24), forming the second bridge (100–10−1 min). It should be noted that there exists another route to give an incomplete cage (IC) (2,3,2) in small numbers (#25), the presence of which was not emphasized as an intermediate in the original nk analysis4 (Fig. 4(b–d)).

• The final stage of the self-assembly shows sequential intramolecular ligand exchanges (∼101 min), reaching the final cage Pd214 (∼101 min) (Fig. 4(d and e)).

The overall self-assembly mechanism is shown in Fig. 1(b), where the major and the minor reaction pathways are indicated with thick black and thin gray arrows, respectively. It is worth noting that non-negligible amounts of a [Pd213X2]2+ IC (2,3,2) (#25) are transiently produced. The experimental nk analysis for the Pd214 cage showed an apparent simple temporal change along the line k = 0.5 (see Fig. 3(c)), so the reaction route through (2,3,2) was not emphasized, simply because k = 2/3 ∼ 0.66 for this transient species. However, as shown in Fig. 3(c) there should be a distinct possibility that even by going through (2,3,2) the nk plots can stay around 〈k〉 = 0.5 because of the average with residual species which have smaller k values. The experimental observation of the signals for (2,3,2) by ESI-TOF mass spectrometry also supports the existence of the IC.4

3.3 Rate-determining step as a result of comparable chemical reactions

In addition to the fact that it was successful to specify the main reaction path from the reactants to the cage with our minimal reaction network, we should note that the numerically obtained rate constants k3 and k4 have comparable values. We interpret this numerical result as follows. It is well known that the ligand exchange reaction of the Pd(II) ion proceeds via an five coordinate transition state, as shown in Fig. 5. This characteristic should necessarily lead to the fact that the rigidity of ditopic ligand 1 significantly decreases the rate of the intramolecular ligand exchanges, due to the restricted geometry of species determined with preexisting bridge(s) between two metal ions. In our calculations, this is clearly reflected in the values of numerically estimated rate constants, i.e., not only k4 but also k3. In the original QASAP study, only the final stages of the self-assembly were identified as the rate-determining steps. However, the present numerical analysis shows that the final step of the reaction network becomes the rate-determining step without considering only the final step to be extremely slow.
image file: c8cp06102e-f5.tif
Fig. 5 (a) Schematic representation of the ligand exchange mechanism on Pd(II) centers, which adopt a square-planer coordination geometry. The associative ligand exchange leads to five-coordinate intermediates and a transition state. (b) An example of the intramolecular ligand exchange reaction, (2,2,5) to (2,2,4). In the formation of (2,2,4) from (2,2,5), the free pyridyl group in (2,2,5) should approach the axial site of the Pd(II) center to form a square pyramidal intermediate and a trigonal-bipyramidal transition state, where the ditopic ligand is highly strained. This causes slow ligand exchanges even though this process is an intramolecular reaction.

4 Discussion

From the above numerical results, it was found that our simple model with only eight rate constant parameters can reasonably reproduce the existence ratios of reaction products and the nk plots obtained experimentally. Additionally, the time evolution of each participant species could be tracked with the stochastic approach applied. Also found was that there can be another route to the precursor species not emphasized by the original QASAP study. Of course there is a possibility that the rate constant of the intermolecular ligand exchange for the numerically specified reaction [Pd213X2]4+ → [Pd214X]4+ + X is slightly overestimated, because our simplified model assumes that all the intermolecular ligand exchange reactions take place with exactly the same rate constant value k1, regardless of which reactions are mainly occurring in the overall reaction network. However, considering that arbitrarily blocking the path from (2,3,2) to (2,4,1) completely leads to an unnatural imbalance of the reaction network (see the network in Fig. 2), the notion seems acceptable that the final cage should partly be formed through the (2,3,2) IC.

Moreover, from our consideration of the intramolecular ligand exchange reaction with rigid ligands, it was revealed that the rate constant of a particular single step in a set of chemical reactions does not need to be extremely small in order for that step to become rate-determining. This conclusion raises another question about the concept of the rate-determining steps in chemical reactions, especially molecular self-assembly possessing a complicated reaction network like the one treated in this study.

As is expected easily, there remain some differences between the results obtained experimentally and numerically, which would be partly due to neglecting larger intermediates with three or more Pd(II) centers. The fact that (3,5,2), (3,6,0), and (3,6,1) signals are found in the ESI-TOF mass spectroscopy4 tells us that it becomes necessary to consider larger species in order to explain the mechanism of this particular cage formation reaction with further details. In principle, it is straightforward to extend the reaction model to include reactions with such large species. However, the reaction network studied here shows not so much evidence of the appearance of large species, so our simple model can be thought of as well reproducing and qualitatively explaining Pd2L4 cage formation. Such systems as (i) other Pd2L4 cages which have been experimentally known to be formed via a route inevitably including larger intermediate species, (ii) self-assembly recognizing chirality, and (iii) coordination self-assembly with more metal centers (Pd3L6, Pd4L8, Pd6L12, etc.) are very important to obtain deeper insight into molecular self-assembly processes. Numerical analyses for those systems are presently underway and will shortly be reported elsewhere.

Conflicts of interest

There are no conflicts to declare.


This work was supported by JSPS Grants-in-Aid for Scientific Research(C) (16K05511), Grants-in-Aid for Scientific Research on Innovative Areas “Dynamical Ordering of Biomolecular Systems for Creation of Integrated Functions” (25102001, 25102002, and 25102005), and The Asahi Glass Foundation. The authors thank Dr Tatsuo Kojima and Mr Tomoki Tateishi for helpful discussion.


  1. S. Hiraoka, Chem. Rec., 2015, 15, 1144 CrossRef CAS PubMed.
  2. S. Hiraoka, Bull. Chem. Soc. Jpn., 2018, 91, 957 CrossRef CAS.
  3. S. Hiraoka, Isr. J. Chem. DOI:10.1002/ijch.201800073.
  4. S. Kai, V. Marti-Centelles, Y. Sakuma, T. Mashiko, T. Kojima, U. Nagashima, M. Tachikawa, P. J. Lusby and S. Hiraoka, Chem. – Eur. J., 2018, 24, 663 CrossRef CAS PubMed.
  5. S. Kai, M. Nakagawa, T. Kojima, X. Li, M. Yamashita, M. Yoshizawa and S. Hiraoka, Chem. – Eur. J., 2018, 24, 3965 CrossRef CAS PubMed.
  6. S. Kai, S. P. Maddala, T. Kojima, S. Akagi, K. Harano, E. Nakamura and S. Hiraoka, Dalton Trans., 2018, 47, 3258 RSC.
  7. M. Nakagawa, S. Kai, T. Kojima and S. Hiraoka, Chem. – Eur. J., 2018, 24, 8804 CrossRef CAS.
  8. T. Tateishi, T. Kojima and S. Hiraoka, Chem. Commun., 2018, 1, 20 CrossRef.
  9. Y. Tsujimoto, T. Kojima and S. Hiraoka, Chem. Sci., 2014, 5, 4167 RSC.
  10. T. Tateishi, T. Kojima and S. Hiraoka, Inorg. Chem., 2018, 57, 2686 CrossRef CAS PubMed.
  11. S. Kai, T. Shigeta, T. Kojima and S. Hiraoka, Chem. – Asian J., 2017, 12, 3203 CrossRef CAS PubMed.
  12. T. Tateishi, S. Kai, Y. Sasaki, T. Kojima, S. Takahashi and S. Hiraoka, Chem. Commun., 2018, 54, 7758 RSC.
  13. T. Tateishi, W. Zhu, L. H. Foianesi-Takeshige, T. Kojima, K. Ogata and S. Hiraoka, Eur. J. Inorg. Chem., 2018, 1192 CrossRef CAS.
  14. A. Baba, T. Kojima and S. Hiraoka, J. Am. Chem. Soc., 2015, 137, 7664 CrossRef CAS PubMed.
  15. A. Baba, T. Kojima and S. Hiraoka, Chem. – Eur. J., 2018, 24, 838 CrossRef CAS.
  16. T. Iioka, S. Takahashi, Y. Yoshida, Y. Matsumura, S. Hiraoka and H. Sato, J. Comput. Chem. DOI:10.1002/jcc.25588.
  17. D. T. Gillespie, J. Comput. Phys., 1976, 22, 403 CrossRef CAS.
  18. D. T. Gillespie, J. Phys. Chem., 1977, 81, 2340 CrossRef CAS.
  19. D. T. Gillespie, Phys. A, 1992, 188, 404 CrossRef CAS.
  20. D. T. Gillespie, Annu. Rev. Phys. Chem., 2007, 58, 35 CrossRef CAS PubMed.
  21. Y. Matsumura, S. Hiraoka and H. Sato, Phys. Chem. Chem. Phys., 2017, 19, 20338 RSC.

This journal is © the Owner Societies 2019