Computational modelling of the steps Involved in photodynamic therapy

Photodynamic therapy (PDT) is a branch of phototherapy that has seen a surge of interest in the last few decades, due to its potential in the treatment of various cancers, infections and heart disease.(Bonnett, 2000) This chapter aims to give an overview of the various photochemical steps involved in PDT as a cancer therapy, and in particular the challenges and insight gained from their theoretical description. After a brief review of PDT in general, in a biological and chemical context, the photochemical steps involved will be discussed, detailing the computational techniques required to model these chemical pathways theoretically. We will detail the methodologies that can currently be applied, as well as their limitations of use at present, and areas requiring further development.


Introduction
Photodynamic therapy (PDT) is a branch of phototherapy that has seen a surge of interest in the last few decades, due to its potential in the treatment of various cancers, infections and heart disease.(Bonnett, 2000) This chapter aims to give an overview of the various photochemical steps involved in PDT as a cancer therapy, and in particular the challenges and insight gained from their theoretical description.After a brief review of PDT in general, in a biological and chemical context, the photochemical steps involved will be discussed, detailing the computational techniques required to model these chemical pathways theoretically.We will detail the methodologies that can currently be applied, as well as their limitations of use at present, and areas requiring further development.

Review of PDT
Photodynamic therapy is perhaps best introduced by considering the photochemical pathways illustrated in the example Jablonski-type diagram shown below, which describes the various processes that would take place within the vicinity of a target cell.
The photochemical processes involved in PDT starts with irradiation (A).The irradiation can take place externally, normally from a laser source, as well as internally, through the use of optical fibers, and targets tissue containing a photosensitiser molecule (B).The sensitiser is functionalised so that it will accumulate in the cells that are targeted, for example in a cancer lesion or an area of diseased skin cells.It is also designed so that it will effectively absorb light of a specific wavelength.The absorption process (C) takes place either through the absorption of one photon per sensitiser molecule, one-photon absorption (OPA), or, if a light source with high fluency is used (e.g., from an intense laser beam), two photons per molecule (TPA).In the Jablonski diagram the eigenstates are separated according to their energy as well as spin multiplicity, where S denotes a singlet state and T a triplet state.These multiplicities, and their relative order, are typical for many organic chromophores that are used as photosensitisers, and the transitions taking place between them are governed by the various section rules for electronic transitions (e.g.non-parity conserving transitions for OPA and parity conserving for TPA).After very fast non-radiative relaxation processes, collectively termed internal conversion (IC) (D), the absorption results in the generation of the sensitiser in its first excited singlet state, S 1 .The sensitiser molecule must Scheme 1. Overview of the photochemical pathways taking place in photodynamic therapy, described through the energy flow starting from irradiation of the photosensitiser molecule and ultimately ending in cell death.then be able to undergo efficient intersystem crossing (ISC) (E) that generates the first excited triplet state of the molecule, T 1 .The sensitiser has to be designed so that ISC takes place with high quantum yield.Any radiative processes, where the molecule will return to its ground state after emitting a photon, are looked upon as loss processes in PDT and should be not effective deactivation channels in an ideal photosensitiser molecule.The same applies to any intra-molecular reactive channels in the singlet or triplet manifolds, which generate new photoproducts.The comparably long-lived triplet state can react in two ways when confined in a biological environment: Type I reactions, which involves the creation of radicals via electron or hydrogen transfer reactions, or Type II reactions, where an energy transfer process (F) takes place from the triplet state of the photosensitiser, and generates the cytotoxic singlet oxygen species (G) from ground state oxygen which is always abundant in a biological environment.During the energy transfer process the photosensitiser is simultaneously brought back to its singlet ground state where it can, in principle, take part in further sensitisation cycles.The singlet oxygen created is exceptionally reactive, and is therefore very short-lived in a biological medium where it has been shown to be toxic to essentially any biomolecule and organelle within the cell.Both pathways ultimately lead to oxidative degradation of the biomolecules within the cell, but in most cases Type II mechanisms are dominant, and will be the pathway considered in this chapter.Only in environments with low oxygen levels are Type I mechanisms observed, which means it will mainly have to be considered for the targeting of hypoxic parts of fast-growing tumours.(Bonnett, 2000;MacDonald & Dougherty, 2001;Missailidis, 2008) The knowledge of the benefits of sunlight on the general health have been traced back over 5000 years with references from Egypt and India as well as ancient Greek texts mentioning Herodetus and Hippocrates recommending "heliotherapy" to their fellow Greeks and documenting the benefits of sunlight on bone growth.(Bonnett, 2000) The development of PDT as we know it today has been in progress since the early 1900s, an example being the 1903 Nobel Prize in Physiology or Medicine awarded to the Danish scientist Niels Rydberg Finsen.(Nobelprize.org, 2011) Finsens work was the first of its kind where light was used actively and directly as a treatment of disease.The first descriptions of the use of a chemical photosensitiser compound in light therapy are also from the early 1900s through the observations of Raab as well as Jasionek and Von Trappen.(Ethirajan, Chen et al., 2011) Raabs work involved the investigation of bacterial cultures and the effect of an added porphyrin dye.He observed completely unexpected cell death in the experiments taking place during the daytime, but none during his evening experiments, which lead him to draw conclusions about the light activation of his dye molecules and their therapeutic outcome.(Allison, Downie et al., 2004) One of the most notable examples from this period is the 1913 experiment by the German scientist Meyer, who injected himself with 200 milligrams of haematoporphyrin, a carbonyl functionalised porphyrin derivative (Fig. 1).After sunlight exposure he suffered extreme swelling and experienced photosensitivity for several months after the experiment.(Dolphin, Sternberg et al., 1998) A more recent discovery in the development of PDT is that a heamatoporphyrin derivative has tumour localising properties as well as being phototoxic.This discovery prompted a spur of development between 1960s and 1980s which ultimately lead to the first clinically approved drug Photofrin®, based on the haematoporphyrin derivative, as well as being the starting point for the development of PDT as a promising tool in cancer treatment.(Hasan, Celli et al., 2010) The current use of PDT in cancer treatment is typically palliative and for very advanced stages of bladder, lung and oesophageal cancers.Its only success as a general treatment is in the case of age-related macular degeneration, which is a leading cause of blindness in the western world, and this still represents the only area where PDT is the first treatment of choice.(Schmidt-Erfurth, 2000) The main limitation in the general application of PDT as an anticancer therapy is the poor level of tissue penetration of the light used.Skin cancers are the only type that can be treated with relative ease, whilst for lung and oesophagus tumours the administration of light normally requires the use of an endoscope.A further problem arising if the light dose can be administrated is the targeting of the photosensitiser compound.Poor targeting leads to general photosensitivity in patients for weeks after the treatment, and even high strength light bulbs can be problematic.General side effects from the treatment of internal cancer are nausea and fever, caused by swelling and inflammation of the nearby tissue.(Davids & Kleemann, 2011) Despite these side effects and even though developments in PDT have not yet resulted in its general use as an anti-cancer therapy, it is still an exceptionally attractive goal and research continues today along many promising avenues.The attractiveness of PDT as a clinical therapy comes from several aspects that could in principle elevate it above current modes of treatment.It is a non-invasive therapy that has no toxicity in the absence of a light source, and compared to surgery it has excellent outcome with respect to scarring and other cosmetic considerations.As the technique generally requires light in the ultraviolet to infrared region the radiation damage is minimal, compared to other forms of radiotherapy, and it has also been suggested that PDT can trigger and help mobilise the immune system against the tumour.(Brackett& Gollnick, 2011) A further major advantage PDT has, in comparison to various forms of chemotherapy, is the inability of the body to develop resistance to the treatment as there is not necessarily a specific target biomolecule or organelle involved in the mechanism of its cytotoxicity.These facts have stimulated continued research as we detail below.One of the most important aspects is a deeper understanding of the spectroscopic and mechanistic features as shown in Fig. 1.First principles computation has only recently matured to the point where it can be routinely applied to some of these problems, and still other challenges remain.Nevertheless such studies are now regarded as invaluable in rationalising the existing classes of photosensitizers, and can now be used in the design aspect of photosensitisers to have a set of desired molecular properties (Arnbjerg, Jimenez-Banzo et al., 2007;Arnbjerg, Paterson et al., 2007;Bergendahl & Paterson, 2011;Johnsen, Paterson et al., 2008).
Any molecule that show promising behaviour for the use as a photosensitiser will have to be further optimized extensively in order to consider targeting within the human body, toxicity, pharmacokinetics, ease of total synthesis of the molecule, and so on.With computational chemistry having reached this stage, systematic studies of molecular structure for PDT properties can hopefully aid in the design of drug molecules for application in this field.

Computational techniques
In order to fully explain the fate of an electronically excited species theoretically the Jablonski diagram above needs to be expanded, and the energy of each state needs to be related to the molecular structure of the system under investigation.This is done by the expression of the energy as a function of the degrees of freedom in the molecule.For any molecular system, apart from simple diatomic molecules, this will result in a multidimensional potential energy surface (PES).This concept, where the motions of the atomic nuclei are treated separately from the motion of the electrons, is the result of one of the fundamental approximations employed in molecular quantum mechanics: the Born-Oppenheimer (BO) Approximation.Conceptually it is based upon the idea that the relatively large velocity and low mass of the electrons in comparison to the nuclei means their motion can be decoupled and treated separately, as the dynamics of these species will take place on vastly differing timescales.The BO approximation can be used with good results in the theoretical description of most areas of chemistry, and the topology of a PES can be used routinely to analyse the thermodynamics and kinetics of most thermal reactions.We will see however that this approximation fails when describing certain photochemical pathways where two PESs can be looked upon as crossing each other, as in IC or ISC when discussing PDT, a problem that needs to be addressed for an accurate theoretical description.Thus, such treatment of so-called vibronically coupled states is one of the most challenging areas, but also one of the most important in Fig. 1 above, as this determines the exchange of the initial photo-energy between the electronic and nuclear degrees of freedom.

Modelling of the electronic ground state
When attempting to describe theoretically the PES of the ground state of a large molecular system from first principles it is convenient to divide the solution to the electronic Schrödinger equation, the total electronic energy, E el , into various part: where E T is the kinetic energy of the electrons, E V is the coulomb interaction between the electrons and the nuclei and the coulomb term E J describes the energy the electrons would have if they moved independently and if each electron also repelled themselves (i.e., the repulsion of the "smeared out" charge clouds of each electron).The final two terms correct for the false parts of the latter, E X , the exchange term, corrects for the fact that electrons with the same spin will avoid each other strongly as a result of the Pauli principle, whilst E C , the correlation term, corrects for the smaller correlation of the movement of electrons with opposite spin.(Gill, 1998) Most computer codes will calculate an approximate solution to these parts of the total electronic energy separately, and the way in which they deal with each of them is what distinguishes various methods from each other.That, and of course the computational cost of each method (in terms of CPU time).

Wavefunction methods
The exact solution to the electronic Schrödinger equation is an intimidating task for any system with more than one electron due to the complexity of dealing with interacting electrons.The first approximation generally applied is therefore to make a crude assumption that the total many-electron wavefunction, Φ, can be constructed as a product of functions, ϕ, which describes one electron at the time, i.e. molecular orbitals.This is the approach used in the Hartree-Fock method, which is the best approximation available when using one-electron molecular orbitals to describe a single electronic configuration, i.e., the orbital occupancy of the system.To account for electron spin and indistinguishability the orbitals are extended to spin-orbitals and described as a Slater determinant, which satisfies the Pauli principle for electrons.The total many-electron Hamiltonian is then used with these one-electron orbitals, and the resulting energy is optimised variationally.Crucially however, the resulting operator after this treatment is expressed in terms of the molecular orbitals obtained as its eigenfunction, creating a non-linear problem that needs to be solved iteratively through a so called self-consistent field (SCF) calculation.(Cramer, 2004;Jensen, 2007;Szabo & Ostlund, 1996) It should be noted that in the HF method the last term in eqn. 1 above is completely unaccounted for, and thus this MO method is generally described as an uncorrelated method.The electronic configuration described by the HF method need not be that of the electronic ground state but if a configuration of a different symmetry is constructed then the variational nature of the method can be used to determine the optimum molecular orbitals for that configuration.Thus, for an excited state that can be qualitatively described by a single electron configuration, of a different symmetry to the ground state, then this method is capable of describing such a state.In this so called Δ-SCF approach two problems arise, (i) most excited electronic states are poorly represented by a single configuration, and (ii) this method is unable to describe excited states of the same symmetry as a lower lying state due to variational collapse.
One key part in carrying out any computation is the choice of basis functions used to describe the one electron orbitals, ϕ, through an expansion of the molecular orbitals via known atomcentred functions (i.e., the well known LCAO expansion).It is a dilemma to choose functions that accurately describe the electron distribution in a given region of space, but at the same time are not so numerous as to make the computational cost too steep.Strategies that are employed to increase the quality of a given basis set are to include more functions of each angular momentum type that are used, as the addition of diffuse functions.(Papajak, Zheng et al., 2011) Diffuse functions have a larger radial component than the valence functions and can therefore more accurately describe wavefunctions at a larger distance from the nucleus, which is useful when considering properties of excited states, as well as properties such as polarisabilities, or when dealing with negatively charged species.A veritable plethora of basis sets exists for excited state computations.
The HF method exactly accounts for the fact that electrons with the same spin will avoid each other strongly, so-called electron exchange, through the anti-symmetry of the Slater determinant.However it completely ignores the smaller correction from the movement of electrons with the same spin, this is termed electron correlation.This missing electron correlation energy is normally only a fraction of the total energy, but it is often crucial for the correct description of a system as it varies greatly with molecular geometry.This is especially true for excited states, where there might be more than one electronic configuration needed to describe the true bonding situation in the system, so called static, or non-dynamic correlation.Thus, static correlation describes the cases where a qualitatively correct wavefunction requires several electronic configurations (i.e., a multiconfigurational wavefunction).This can either be due to spatial quasi-degeneracies (i.e., small HOMO-LUMO gaps) or to spin-symmetry (vide infra).The correlation that takes the instantaneous Coloumb repulsion between electron pairs into account, rather than the mean field SCF orbital description, is termed dynamic correlation.This can either be recovered from single or multiconfigurational reference states.
The problem of accounting for the missing electron correlation energy is treated differently in various methods, but the common aspect of all approaches is the understanding that electron-electron interactions can be described more accurately if more than one configuration is used to express the final many electron wavefunction, Φ. (Jensen, 2007) Conceptually the most straightforward solution to the correlation problem is to generate more electronic configurations by promoting electrons in occupied orbitals to virtual, unoccupied, orbitals, so that the system is described by a linear combination of all excited configurations.This is termed the Configuration Interaction (CI) method: 0 ...
If the optimisation of the CI expansion coefficients (also known as amplitudes), a i , were to be carried out over all possible configurations the CI wavefunction would be the exact wavefunction, and hence give an exact solution to the time-independent Schrödinger equation within that given one-electron basis set.This is an extraordinarily difficult task as the number of possible configuration will grows factorially with system size, and the full CI wavefunction normally needs to be truncated in order for it to be routinely used.Even though a truncated CI wavefunction will accurately account for a large fraction of the missing correlation energy, there is one unfortunate drawback.Truncated CI expansions are not size extensive, which means that as the system size gets larger a proportionally smaller amount of the correlation energy is returned at a given truncation level.A related issue is that the CI method is also not size consistent, i.e. the energy is not proportional to the number of non-interacting units comprising the system.This is a significant flaw when calculating for example dissociation energies, as the energy of two dissociated molecular fragments will not be equal to the energy of the two fragments calculated separately, as some multiple excitations are omitted.(Piela, 2007) A method related to CI is the Coupled Cluster (CC) method.It is arguably the most sophisticated and accurate method for retrieving some of the dynamic correlation energy in use today.The cluster operator T is related to the exact wave function through a reference determinant, usually the HF determinant, through the coupled cluster ansatz; where the CC wavefunction is generated through an exponential expansion of the cluster operator that generates the substituted configurations from ϕ in the same way as in CI theory.
It is assumed that the cluster operator is a sum of excitations operators for all possible excitations from occupied to unoccupied orbitals; where T1 is an operator for single excitations and so on.If all possible excitations are considered the result will be the same as for a full CI treatment, and again this full treatment is too computationally expensive to be practical.Yet again a truncation of the wavefunction must take place, and this is where the difference to CI becomes evident.Unlike CI methods the excitation operators will automatically include important so-called disconnected terms, so truncating at T2 in the cluster operator automatically includes T2 2 , T2 3 up to T2 N , from the Taylor expansion of e T , which results in not only a more accurate description of electron correlation but also generates wavefunctions that are both size extensive and size consistent.(Piela, 2007) The disconnected terms like T2 2 can be understood to account for two simultaneous pair-wise interactions away from the HF-SCF mean-field description.The one disadvantage is that the CC method is not variational (the exponential ansatz unfortunately does not generate truncated left and right states in an expectation value of the energy).(Crawford & Schaefer, 2000) However this is not as important as the advantages gained and therefore CC theory occupies the prime position in correlated wavefunction theory.A final advantage the CC approach shares with CI theory is that as the excitation operator is expanded one can systematically improve the properties of the molecule calculated from the resulting wavefunction (with these methods for ground state properties, but see below for excited state properties).However, even with truncation, the computational cost of CC computations are very high and the size of the molecules that can be investigated are limited, and specialist techniques, such as the use of significant computer resources as well as the development of non-iterative versions of the methods, need to be employed if larger systems are to be investigated.(Kowalski, Krishnamoorthy et al., 2010) In order to routinely be able to probe photosensitiser molecules the size of substituted porphyrinic macrocyclic systems different approaches needs also to be considered, for example through the use of Density Functional Theory (DFT), although for the core chromophores of various photosensitisers coupled cluster theory has set important benchmarks against which lower cost methods can be calibrated.

Density based methods
In DFT an alternative starting point for the calculation of the total electronic energy is employed through the use of the total electron density, ρ.The main idea of DFT is that the total energy can be completely determined by the total electron density of the system through a functional that connects it with the energy.(Koch & Holthausen, 2001) The advantage to using the density, rather than the molecular wavefunction, is that the problem will not depend on the number of electrons present in the system, and therefore avoid becoming increasingly complicated as the system increases in size (i.e., the energy ultimately depends on only 3 spatial coordinates rather than 3n for an n-electron wavefunction).In order to develop both an accurate and computationally tractable scheme W. Kohn and L.J. Sham developed a fictitious model system where the electrons do not interact at all so that they are described exactly by a set of molecular orbitals.Crucially, they realised that the density for this non-interacting system can be exactly the same as that of the same system with interacting electrons, and this was therefore used to develop the so-called Kohn-Sham expression for the ground state electronic energy (where each term is obtained from an expression involving molecular orbitals): As the expression deals with non-interacting electrons, the kinetic energy term is not complete since the exact electronic kinetic energy cannot be written as a functional of the density, and its correction is introduced in addition to the exchange and correlation terms, which are all then grouped together in the E XC KS term.All other terms can be evaluated exactly, and so the main problem of DFT is reduced to the estimation of the exchange correlation energy.The way in which various DFT functionals retrieve the exchange correlation energy varies, and is what sets them apart.A number of strategies exist and the improvement of various functionals is generally derived from both theoretical foundations by expanding the parameter set linking the energy with the density (via gradients etc.), and the fitting of such parameters to experimental data.
The Kohn-Sham equation is derived from the Shrödinger equation for the non-interacting system and its solution yields the so-called Kohn-Sham orbitals, φi (r), which are defined such that: These orbitals are then expanded in a basis set, and arranged in a single Slater determinant, both of which are completely analogous to HF MO theory.When finding the energy of a specific molecular structure using DFT the Kohn-Sham equation is solved self consistently using an electron density that initially is a guess in order to calculate E XC KS , and hence the energy of the optimised system.(Koch & Holthausen, 2001)

Multiconfiguration methods
The problem with correlation that we have dealt with so far takes the instantaneous Coloumb repulsion between opposite spin electrons into account, the dynamic correlation (since the Kohn-Sham non-interacting system has a single configuration wavefunction).For some systems however there is also a further complication in that more than one electron configuration might be needed to fully describe the bonding situation in the system.This is especially true for most open-shell systems, of which the oxygen molecule is a prime example, as well as in the description of almost all excited states of closed-shell systems.This effect is the static or non-dynamic correlation discussed above, and is clearly something we need to be able to model if we are to accurately describe all photochemical pathways involved in PDT.In the cases where static correlation is important the use of one single determinant to describe a many electron wavefunction, with the correct spin, is not always possible.There has however been some single-reference methods developed recently that can describe some multiconfigurational situations, for example the so-called spin-flip approach.In this approach a high-spin triplet reference state is chosen, M s = 1.The target state, e.g., the M s = 0 state of an open-shell singlet state, can then be described as a spinflipping excitation from the reference state.(Shao, Head-Gordon et al., 2003) This is a promising avenue of attack for the application of DFT to open-shell systems.A more common approach in order to deal with static correlation however, is to turn to a strategy of taking a linear combination of a number of electron configurations (as in CI theory described above), and define a multiconfigurational wavefunction.This wavefunction is then optimised variationally by an SCF procedure in which all expansion coefficients are variables.That means that the expansion coefficients for the basis functions (the LCAO parameters), which describe the molecular orbitals, as well as the configurational expansion coefficients will be optimised, which adds flexibility to the wavefunction.Thus the optimum orbitals for the given CI-type wavefunctions are obtained, unlike in standard CI theory where the reference orbitals are fixed.The most popular variant of MCSCF theory is the complete active space (CAS) SCF approach.The popularity stems mainly from the fact that the choice of which electron configurations that are important to include in the MC wavefunction can be difficult, and this is automatically taken care of in CASSCF.In the CASSCF approach an active space is selected, that includes specific orbitals, important to describe the chemistry one is interested in, e.g., for benzene photochemistry an obvious choice is the six π-electrons and six π-orbitals.(Li, Mendive-Tapia et al., 2010) In aromatic organic molecules in general, with or without heteroatoms, the frontier molecular orbitals, FMOs, (the Highest Occupied Molecular Orbitals, HOMOs, and the Lowest Unoccupied Molecular Orbitals, LUMOs) are crucial to describe any photochemical mechanisms, and needs therefore to be described accurately.Orbitals of less importance, that are expected to stay either fully occupied or un-occupied in all determinants, are kept unchanged, and thus no configurations arising from changes to these doubly occupied or empty orbitals are included.The remaining orbitals are considered to be active, and a full CI is carried out on this subset.CASSCF is exceptionally flexible and can be used to compute essentially any part of the PES of excited as well as ground states, including surface crossings, which are described later.
Fig. 2. Overview of the CASSCF strategy, where a full CI is carried out in the active space, whilst the occupied and unnoccupied orbitals are kept constant in all configurations.

Modelling of excited states
The single reference methods mentioned so far are all constructed for the description of the ground state potential energy surface of a molecular system, in particular around regions describing stable configurations of the molecule (i.e., PES minima).When dealing with photochemical events obviously more than one PES needs to be considered.
There are many approaches to the theoretical modelling of the properties of a molecular system, such as excitation energies, one example being state specific methods.In these a trial wavefunction is chosen so that it corresponds to the state we are interested in, which could be an excited state rather than the ground state.(Nicolaides, 2011;2011) This is the idea behind the Δ-SCF method discussed above.More powerful are methods that can construct arbitrary excited state wavefunctions of any symmetry.These are most often obtained as higher roots in CI-type eigenvalue problems, and as such CI methods are the method of choice in constructing both qualitative and quantitative excited state wavefunctions.For example the simplest qualitative correct excited state wavefunctions can often be represented by the CIS approach.Also in this manner CASSCF can generate excited wavefunctions of the same accuracy as for ground electronic states.It should be noted that such excited states invariably always require such a multiconfigurational wavefunction at most points on a photochemical reaction path.(Malmqvist and Roos, 1989;Olsen, 2011) Dynamic correlation effects can then be added to such multiconfigurational wavefunctions by using the underlying multi-reference configurations for further CI (MR-CI), or perturbation theory treatments (e.g., CASPT2).(Pulay,2011)

One-and two-photon absorption
One method for the determination of excited state information that has been very successful and popular in recent years is the so-called response function method.The reason for its popularity is that it can be adapted to suit DFT, MCSCF as well as CC methods.Thus, most molecular properties, such as excitation energies and transition moments can be obtained from a response function calculated from a ground state quantity (the wavefunction in MCSCF and CC theory, and the density in DFT).In response theory any molecular property is expressed in terms of a response function, as a property of a molecule can be defined as the response of an expectation value of an operator to a perturbation.(Bast, Ekstrom et al., 2011) When dealing with perturbations in the form of an external field interacting with a molecular system, as when dealing with absorption of light, a time-dependent approach needs to be considered.The time-dependent perturbation term is expressed by an interaction operator, which is added to the Hamiltonian for the problem.
The linear response function for operator P with respect to perturbation Q can be written as: An example is the linear response of a system subject to a homogenous electric field of frequency ω b , where the perturbation operator then will refer to the electric dipole operator, μ .This equation is identical to the dipole-dipole polarizability tensor at frequency ω b , i.e. the frequency-dependent polarizability.This has poles (i.e.where the denominator is equal to zero) at the excitation energies of the system.The residues of the response function (i.e. the numerator at the poles) give the transition moment (related to the probability of a transition between states).The poles of the response functions are obtained as solutions to an eigenvalue equation, with the eigenvectors giving the contribution from orbital excitations to each excited state.This method is very efficient since there is no need to construct all the excited states wavefunctions explicitly as in the sum-over-states expression in the equation above.All the information regarding the excited state is contained in the response function, so there is only a need for one well defined wavefunction or density, such as a ground state wavefunction from DFT, MCSCF or CC, or density from the Kohn-Sham orbitals in DFT, in order to access information on the excited states.Thus the response function approach allows us to access excited state information from nominally ground state methods and take advantage of what they bring (i.e., systematic accuracy in CC theory, low cost in DFT) As previously mentioned, there is also a possibility of a two-photon absorption process for the excitation of the photosensitiser, as per Scheme 2. This has numerous benefits for the application to PDT, such as access to the so-called tissue transparency window (i.e., avoiding absorbance from water, melanin and haemoglobin) and increased spatial resolution for the treatment, thus the accurate modelling of this non-linear absorption process is of great importance in the development of ideal photosensitiser molecules.
Scheme 2. Overview of the TPA process, normally described as taking place through a virtual state, v.
The use of a state specific method for the computation of TPA is inaccurate, as a few state model would need to be constructed, so the general strategy is the use of response theory through higher order terms of the response function.For TPA the quadratic response function needs to be considered under two dipole perturbations, one of frequency ω b and one of frequency ω c .The energy representation of the generalised quadratic response function is then described as: When the operators P, Q and R are dipole operators then this expression is identical to the first dipole hyperpolarisability tensor at frequencies ω b and ω c .This frequency-dependent hyperpolarisability can describe non-linear phenomena involving two photons of frequencies ω b and ω c respectively, as well as TPA information, where two identical photons are absorbed simultaneously.(Paterson, Christiansen et al., 2006) The residues of the quadratic response function can be shown to give the conventional expression for the twophoton absorption transition matrix element, M, between states i and j through a virtual state (written as linear combination of all real eigenstates v): where ω the frequency of the irradiating light, and ω v is the excitation frequency of state v.
The labels α and β refers to the Cartesian coordinates x, y and z of the TPA tensor.The transition matrix element is then used to calculate the rotationally averaged two-photon absorption cross section, δ.
Even though the formal energy expressions for the quadratic response function appears to involve a tedious sum over all states of the molecule to describe the virtual state v above, as mentioned previously the actual computational machinery of response theory calculates the properties directly without constructing the states explicitly which saves both time and computational effort.Additionally, one does not need to make any a priori estimate of the importance of certain states to the TPA tensor.Not only that, this technique has also been proven to be exceptionally accurate in the prediction of TPA spectra.When response theory is used with a CC ground state wavefunction quantitative agreement with experimental results can be achieved, and for a DFT ground state density very good agreement is also often obtained.(Arnbjerg, Jimenez-Banzo et al., 2007;Hattig, Christiansen et al., 1998) When it comes to DFT however, the choice of functional is crucial for the application of quadratic response theory, as it has to be able to model excitations to and from the virtual state, which essentially means it has to be considerably more robust than the functionals routinely used together with linear response theory for one-photon absorption.If this requirement is satisfied quadratic response DFT has been shown to be very accurate, and even comparable to high-level CC calculations for TPA properties.(Arnbjerg, Paterson et al., 2007;Johnsen, Paterson et al., 2008;Paterson, M. J., Christiansen et al., 2006) This is very encouraging, considering the size of the photosensitisers we are interested in modelling for the description of TPA in PDT (vide infra).The improvement of non-linear absorbing qualities of photosensitiser molecules can involve many complicated synthetic routes (Drobizhev, Karotki et al., 2002;Drobizhev, Meng et al., 2006) As the benefits of the use of TPA in PDT are plentiful, a theoretical input will be invaluable, as it can give important information on molecular structure-TPA relationships which are crucial for the development of TPA PDT on a daily basis.This has for example been demonstrated in an experimental and theoretical study, where the porphycene macrocycle was shown to be a promising photosensitiser system for TPA PDT.(Arnbjerg, Jimenez-Banzo et al., 2007) We have carried out further computational TPA studies of this porphycene molecule and shown that the non-linear absorption can be effected to a somewhat unexpected degree by the substitution of heteroatoms in the macrocyclic core.(Bergendahl & Paterson, 2011) The structural changes afforded very effective absorption in the so-called Soret-region of the TPA spectra (an unusual feature for the porphyrin photosensitisers so far investigated for TPA PDT) without having a large impact on the one-photon absorption spectra.We investigated a set of dioxaporphycenes, where the pyrrolic nitrogen atoms were substituted for oxygen atoms systematically as per Scheme 3.
After optimising the geometries of the molecules, the OPA and TPA was calculated using linear and quadratic DFT response theory respectively, with the CAM-B3LYP functional and Scheme 3. The optimised geometries of the three dioxaporphycenes investigated, with resulting symmetries.N atoms in blue, and O atoms in red.
6-31G* Pople-type basis set.Previous work has shown that the CAM-B3LYP functional is very accurate at modelling TPA, as it has correct asymptotic behaviour in the exchangecorrelation functional, which gives a good description of absorption to and from the virtual state.(Paterson, Christiansen et al., 2006;Rudberg, Salek et al., 2005) Fig. 3. Summary of the OPA spectra for the dioxaporphycenes investigated.Note the absorption band at 300 -350 nm, the Soret-band, and at 520 -580 nm, the Q-band, which are characteristic bands for the porphyrin family The OPA results are summarised in the simulated absorption spectra in Figure 3, and show clear Soret-(~ 300-350 nm) and Q-bands (~ 500 -600 nm) for all species.As expected, there are no major effects observed on the spectral profile upon substitution of the core heteroatoms, compared with porphycene.The TPA spectra are presented in Figure 4.It is important to note that the absorption wavelength indicated is the total for the two-photon process, i.e. the actual excitation takes place at twice the wavelength.The spectra all show a main absorption in the Soret-region (~ 250 -270 nm), and there is a significant disparity between the calculated TPA cross-sections for the four isomers.These two facts can be explained in terms of one-photon resonance enhancement of the absorption.This phenomenon is best discussed in terms of the sum-over-states expression in the equation for the transition matrix element M, as seen previously, where the magnitude of M depends on the difference between the excitation frequency of the v-th state and the frequency of the total excitation, from the magnitude of the denominator in the expression.As the virtual state can be described as a linear combination of all real eigenstates in the system, a TPA via a virtual state that can be described mainly in terms of a real state will be significantly enhanced.Thus, an allowed OPA Q-band transition accidentally degenerate at half the Soret excitation energy will greatly enhance the Soret TPA transition.These results suggest that the tuning of TPA does not have to be carried out through large structural changes, and small electronic factors can be used to improve the magnitude of the TPA cross-section, hence broadening the scope for the development of TPA photosensitisers in PDT.

Photochemistry on excited state potential surfaces
Even though the information we gain from response theory can accurately predict the electronic spectra for both one and two photon absorption processes, it is insufficient to give a full overview of the whole reaction path.This is true when we are modelling reactive photochemistry in general, and the paths involved in PDT in particular where IC and ISC plays a significant role in the understanding of the reaction path.The modelling of reaction paths which cross multiple surfaces have been studied both experimentally and theoretically for number of years, and the understanding of the pathways involved in the photochemistry of small molecules have expanded in great detail.(Zewail, 2000;2000) The non-radiative IC pathways have been known to exist for a very long time, an early description being Kashas rule, which states that radiative transitions and ISC will take place from the lowest energy excited level of that spin multiplicity due to the very short time-scale of the IC processes between upper electronic levels.(Kasha, 1950) More recently however, the description and understanding of these events have experienced a dramatic increase through the development of electronic structure methods to theoretically model these events.It is known that IC events are normally dominated by so called conical intersections, where two, or more, surfaces cross each other.Conical intersections are mathematical objects that have been known since the 1930s, and many excellent reviews exist on the topic.(Matsika & Krause, 2011;Worth & Cederbaum, 2004;Yarkony, 1998) They were well known early on to exist in high symmetry species, and manifest through phenomena like the Jahn-Teller effect, showing that molecules in orbitally degenerate states undergo non-totally symmetric distortions to remove the degeneracy.(Paterson, Bearpark et al., 2005) However, advances in electronic structure theory in the last 20 years or so have highlighted that accidental conical intersections, not requiring high symmetry, are exceptionally common in polyatomic systems, and especially in heteroaromatic molecules, which is of interest from a photosensitiser point of view.
Fig. 5.When the energies of the crossing surfaces are plotted agains two coordinates making up the so-called branching space, the surfaces meet in a characteristic double cone shape forming a funnel along the reaction path.
The experimental evidence for the presence of conical intersections often require advanced set-ups and generate complex data sets, such that a theoretical input is often required for full characterisation.(Hadden, Wells et al., 2011;Oliver, King et al., 2011) In non-radiative processes, where a system moves from potential energy surface to another, this population transfer is most efficient at points of degeneracy (i.e., conical intersection seams).Furthermore, the timescale for such transitions is frequently found to be ultrafast (i.e., femtosecond) as a single molecular vibration can cause a change in electronic state.Clearly non-adiabatic terms are very important, as the dynamical processes involve the coupling of nuclear and electronic motion to be taken into account if the full reaction path is to be understood.(Worth & Cederbaum, 2004) Modern computational codes can straightforwardly locate stationary points on a surface, but in the case of two surfaces crossing, the actual crossing point is not a stationary point of either of the surfaces.The intersection seams are of dimension 3N-8 for an N-atom molecule (every point on the seam being a conical intersection geometry).Modern optimisation algorithms optimise to points of zero gradient in this so-called intersection space, generating minimum energy crossing points.(Paterson, Bearpark et al., 2004;Yarkony, 2005) These act as photochemical funnels for the reaction in the upper electronic state.It is important to note however that this minimum energy crossing point might not be the only one of importance when it comes to the generalisation of various reaction paths, but a good characterisation of the full seem is crucial.(Paterson, Robb et al., 2005) The procedure that is needed to find the minimum energy crossing point is actually fairly straight forward, at least in principle, and involves computing both the energy and its gradient on both surfaces, as well as the non-adiabatic coupling, which essentially describes the extent to which the electronic states are coupled by nuclear motion (i.e., a break-down of the BO approximation discussed earlier).Multiconfigurational methods are the most common, and appropriate, for the calculations of the IC pathways, as they can handle the multi state character needed near a crossing point.However, as the ubiquitous nature of conical intersections, and their importance in the description of photochemistry, in molecules with biological importance, such as photosensitiser molecules, becomes more evident there is a need to develop methods that can handle systems of a larger size.The use of TD-DFT, which is successful for excited states in the vertical region, in the characterisation of crossing points have been investigated, and has been shown to be of some use for the location of minimum energy crossing points between excited levels (i.e., not with the ground electronic state) They can not yet describe excited state surfaces accurately in general, so care must be taken in their application to problems in reactive photochemistry as opposed to spectroscopy.(Dreuw & Head-Gordon, 2005) When it comes to the description of surface crossings of larger molecules, alternative multiconfigurational approaches are required.This can be carried out in two general ways: the truncation of the computational method we want to use, or the truncation of the molecule to a representative model system.Bearpark et al, investigated a collection of photochemical reactions involving multi-state processes using a CASSCF starting point.(Bearpark, Ogliaro et al., 2007) They primarily considered so called hybrid methods where the molecular system is partitioned so that one area, such as the reaction centre or chromophore, is treated with a more accurate, and computationally more expensive, model.
In addition they investigated a conventional method, where the whole system is treated at the same level of theory, although biased towards a particular region of the molecule.At the moment the practical limit of the use of CASSCF lies in the region of 14 active electrons in 14 active orbitals.Methods that can combine a quantum mechanical method (QM) with a Molecular Mechanics method (MM) (which determines molecular properties through the use of empirically built force fields) are generally referred to as hybrid methods.Hybrid methods, such as ONIOM investigated in this study, rely on the fact that some parts of the molecule are more crucial in the representation of the chemical reaction we want to study, whilst the remainder are viewed as playing a minor role.The ONIOM version of QM/MM method has the advantage of being very general, as in principle more layers of levels of theory can be applied.This method was found to be appropriate for excited state potential energy surfaces, provided the chromophore could be accurately separated from the full molecule.This highlights the importance of these models in the investigations of chromophores in for example a protein or solvent environment.When it comes to molecules where the chromophore is essentially the entire system however, these methods were not always appropriate.This is a key factor when dealing with photosensitisers for use in PDT, as they normally are of a large structure incorporating many aromatic features, such as porphyrins, where the whole system is crucial for the modelling of its photochemistry.For molecules such as this the strategy would be to turn to models where the method itself, rather than the molecular structure, has been truncated, such as the RASSCF approach.(Merchan, Orti et al., 1994;1994) The RASSCF method is analogous to CASSCF where the active space has been divided into smaller parts, as to reduce the total number of electron configurations generated.The active space is separated according to maximum and minimum occupancy criteria, where the RAS 1 space only allows a limited number of vacancies, holes, and the RAS 3 space only allows a limited number of electrons.In the remaining RAS 2 space a full CI is carried out, as in the general CASSCF approach.In this study it was found that the usefulness of the various methods studied was very systemdependent.However, providing the partitioning of the system, or selections of orbitals in the RASSCF case, were adequate the calculations lead to small errors when modelling the excited state surfaces, when compared to CASSCF, on a 14-electron system.(Bearpark, Ogliaro et al., 2007) We are currently investigating the use of RASSCF for analysis of surface crossings in molecules of the size of common photosensitiser molecules.The RAS space chosen here highlights the difficulty in the use of multiconfigurational methods for large molecules, as essentially the whole molecule forms a crucial part of the chromophore (Fig 7).The selection of the orbitals to be included in RAS-2 has to be founded on photochemical background of the problem.In the case of porphyrin derivatives it is well known that the two dominating absorption bands in the spectrum (the Soret-and Q-bands) arise from excitation from the HOMO and the HOMO -1 to the LUMO and the LUMO+1.The use of these FMOs to describe the absorption characteristics of porphyrin systems is termed the Gouterman Four Orbital (GFO) model and it is obvious that an accurate description (i.e., full CI in the RAS 2 space) of these orbitals will be crucial for an accurate RASSCF calculation (Fig 7 ).A methodical description of a 26-electron 24-orbital system like these is already a challenging task, and in order for these methods to be useful for the investigations into photosensitiser molecules with various specialist properties, like large TPA cross section, advances in both electronic structure codes and computational performance will be crucial and lead to important insight for these types of investigations.where the whole π-system have to be accounted for.RAS 2 includes the frontier molecular orbitals, FMOs, of the porphyrin macrocycle, as these have been shown to be vital for an accurate description of the absorption in molecules of this type.
In the computation of the crossing of states with different spin symmetry, ISC, it is necessary to consider spin-orbit coupling, the coupling of the net spin and orbital angular momenta of the electrons.This can be viewed as a small amount of spin being mixed between two states, which lead to intensity from 'spin-allowed' transitions getting transferred to formally 'spinforbidden' transitions, hence facilitating the increased probability of the latter.The process of ISC can be modelled through multiconfigurational methods, as described above for an IC process.It should be noted that for weak spin-orbit coupling in photosensitisers the geometries at which the singlet and triplet surfaces cross are the dominant features that need to be understood.These crossings are 3N-7 dimensional seams in contrast to conical intersections due to the non-adiabatic coupling being zero by symmetry.Standard conical intersection algorithms can be routinely applied to find minimum energy points on these ISC seams.The use of RASSCF to model this is now possible given recent advances in the RASSCF method, in particular the use of gradient driven RASSCF.
There has also been considerable development in the use of DFT to model a change in spin in the ground state of a molecule, mainly spurred from the interest in the modelling of chemical reactions involving a change of spin in transition metal compounds, which can be of interest from the point of view of photosensitiser molecules.(Poli & Harvey, 2003) Even though DFT has had most success as a ground state method, the modelling of crossings between spin-states can be carried out effectively, as the event can be viewed as a single surface reaction.This is done by considering the reactants and products of the photochemical reaction as being located on two diabatic surfaces that have significant mixing in the area of the transition state.If the surfaces of the two spin-states are described accurately by DFT, the method can be used with good results to predict the energy of the crossing point between them.Of course this only applies to the crossings between the lowest energy states of each multiplicity.The main issue with DFT however is still the choice of functional to apply to the problem, as there is no direct recipe even for what class of functional will preform better in these cases.(Reiher, Salomon et al., 2001) Fig. 8. Two diabatic surfaces, each desribing a spin state, crossing in an area which can be viewed as the adiabatic TS.

Photosensitisers in PDT
As in the development of any drug molecule there are many aspects to consider in the advancement of an ultimate photosensitiser molecule, such as pharmacokinetics, targeting and cytotoxicity.However, in the case of a photosensitiser molecule, there is also the absorption qualities to consider and, in particular, the location of the wavelength of maximum absorption in relation to the of tissue transparency window. .(Peng & Moan, 2003) In the body, the most recognisable porphyrin derivative that exists naturally is coordinated to Fe in the haem structure, present in the haemoglobin protein.The unchelated porphyrin structure is known as protoporphyrin IX (PpIX), which has been shown to be an effective photosensitiser.(Kennedy, Pottier et al., 1990) The strategy for PpIX treatment is to administer a synthetic version of δ-Aminolaevulic acid, the biosynthetic precursor to PpIX, which leads to a build-up of PpIX in the targeted area.In Europe, PpIX is approved for use in PDT of actinic keratosis and basal cell carcinoma.(Morton, McKenna et al., 2008) The most common photosensitiser in use today, however, is still the heamatoporphyrin derivative (HpD), which consists of a mixture of oligomeric porphyrin chains.In its purified form it is known as Photofrin®, which was the first PDT sensitiser to be approved for clinical use.(Elsaie, Choudhary et al., 2009) Photofrin or HpD is approved in Europe, the USA, Canada and Japan and is used mainly for cervical, bladder and gastric cancers.
Although Photofrin®, and other porphyrin derivatives, show excellent singlet oxygen quantum yields the longest wavelength at which they are excited is at 630 nm.(Dolphin, Sternberg et al., 1998) This maximum absorption is within the transparency window of human tissue as desired but not at an optimum position.Light of a longer wavelength would penetrate deeper and the clinical use of the photosensitiser would be expanded.This has lead to extended research in ways to manipulate the wavelength of maximum absorption in porphyrin related molecules and has lead to substituted porphyrins being the most studied class of sensitiser molecules.(Peng & Moan, 2003) One strategy to alter the porphyrin structure is to add substituents in the meso-position.Meso-substituted porphyrins were investigated by Oliviera and co-workers, and they could demonstrate that unsymmetrical substitutions on the porphyrin molecule had neglible effect on the absorbance as well as on the fluorescence spectrum.The effect on the singlet oxygen quantum yield was also neglible.These were significant results, as they points out a porphyrinic site that is available for potential modification with drug targeting in mind, without necessarily affecting the chromophore itself.(Oliveira, Licsandru et al., 2009) However, most known porphyrins do not absorb at longer wavelengths than Photofrin, so it is likely that interest as well as importance of pure porphyrins for use a photosensitiser agents in PDT will decrease in the future, whilst concentration is being shifted towards porphyrin derivatives, so called second generation photosensitisers.
Substituting one N atom in the porphyrin core for a heavier S or Se atom, in order to facilitate larger spin-orbit coupling, have little effect on the absorption spectrum, apart from a shift in the absorption maximum to slightly longer wavelengths (~665 nm), (Rath, Sankar et al., 2005) (Detty, Hilmey et al., 2002) however without any effect in the generation of singlet oxygen.Substitution of a further N atom with S or Se moves the absorption maximum to longer wavelengths, but shows a negative effect in the generation of singlet oxygen.(Detty, Gibson et al., 2004) Even though the singlet oxygen generation is decreased in the structures where two core heteroatoms have been substituted, they have been proven to be more effective than meso-substituted pure porphyrin cores in in vivo studies.(Detty, Gibson et al., 2004) The avenue for the use of core-substituted macrocycles in TPA PDT is still in its early stages of investigation, but does show some promise.(Bergendahl & Paterson, 2011) Porphyrin isomers that has one of their β-carbons and pyrrolic nitrogen atoms switched are termed confused porphyrins and they show a maximum absorption at longer wavelengths than porphyrin (~ 730 nm).(Furuta, Asano et al., 1994) Even though these types of compounds have potential as photosensitisers, with moderate singlet oxygen yield, their difficult total synthesis is a major obstacle for their development as photosensitiser agents in PDT.
A molecular structure related to the porphyrin core can be found in the porphycene compounds, a further reduced porphyrin isomer.The porphycene core structure have been shown to experience a large absorption at slightly longer wavelengths than porphyrins (∼ 633 nm) which is thought to be due to lower molecular symmetry, and in vitro studies prove a modest singlet oxygen quantum yield.(Stockert, Canete et al., 2007) Porphycenes have also been shown in in vivo studies to exhibit selective uptake in membrane related organelles, which was proven to ultimately lead to apoptosis of the cells after application of radiation.(Stockert, Canete et al., 2007).The importance of the porphycene system has been suggested to lie mainly in the application of TPA PDT, as early theoretical studies suggests its superiority at effective TPA compared to many related porphyrin derivatives.
As the main focus of the investigations into effective photosensitiser system had been dominated by synthetic developments of porphyrin related systems, this is also the avenue that recent theoretical efforts have gone down.Hence the theoretical methods mentioned in this chapter have been discussed with systems like these in mind.It is however important to keep in mind that non-porphyrin photosensitisers also are being investigated, including various dye molecules and systems on the nano-scale such as semiconductor quantum dots and fullerenes.(Bakalova, Ohba et al., 2004;Mroz, Tegos et al., 2007) As the field of theoretical chemistry develops, with improvements in algorithms as well as advances in computational hardware, the input into PDT research has potential of becoming very important in the future.

Oxygen in PDT
The species we have so far referred to as 'singlet oxygen' in our discussion on PDT is the first excited state of molecular oxygen (a 1 Δ g ).Reactions involving interactions from this state play a key role in many photodegradation processes in photochemistry as well as photobiology and various applications exists, aside from PDT, examples including mechanistic cell death studies, bleaching processes and general organic synthesis applications.(Schweitzer & Schmidt, 2003) Fig. 10.Molecular orbital diagram for oxygen.The occupation of the antibonding π-orbitals, which give rise to the various excited states of O 2, is indicated.Adapted from (Paterson, M. J., Christiansen et al., 2006) The a 1 Δ g  X 3  g  transition is arguably the most forbidden electronic transition in nature, with violation of not only the spin selection rule, but parity and orbital selection rules as well, and in fact all the radiative processes of the singlet states of the O 2 molecule are forbidden to a various degree.This leads to an extraordinary metastability of these states in the gas phase, with lifetimes of 45 min for the 1 Δ g state and 7.1 sec for the 1  g  state.(Kasha & Brabham, 1979) Crucially, however, the degree that these radiative transitions are forbidden decreases rapidly as soon the molecule is perturbed by other atoms or molecules.The nature of this perturbation is generally thought of as a weakly bound collision complex, M-O 2 , where M could be a number of species, including solvents or O 2 itself.In the case of sensitisation O 2 ( X 3  g  ) forms a collision complex with an excited state of a photosensitiser molecule, as will be described more later.The weakness of this interaction presents problems for computational modelling.Firstly, many possible conformations may be possible and some dynamical sampling will be required.There are at present no good forcefields to model this with MD.Secondly, methods such as DFT do not generally describe such non-covalent interactions very well, although empirically corrected functionals describing such dispersion effects have recently been developed.
As illustrated in the molecular orbital diagram in Figure 10, 14 electrons (including the 1s 2 pairs) of the O 2 molecule can be paired up in the lowest lying molecular orbitals.This leaves two remaining electrons that have to occupy a doubly degenerate set of orbitals ( g ).There are several possible ways to accomplish this.This is a basic and simple view of the electronic structure of O 2 , but it highlights some important aspects of the system, such as the fact that the differences in the three lowest lying states in O 2 are only a result of the pattern in which the final two electrons occupy two orthogonal and degenerate molecular orbitals.It also highlights the fact that the photochemically important a 1  g and X 3  g  states are both openshell states, i..e, states involving partially occupied orbitals.This open-shell character is crucial to consider if we are to attempt to accurately model O 2 , and M-O 2 collision complexes computationally.(Schmidt & Schweitzer, 2003) (Paterson, Christiansen et al., 2006) When modelling the oxygen molecule, and its role in PDT, the important states are in fact all of open-shell character to some extent, and often require several electronic configurations to describe them.This open-shell character gives rise to physical phenomena that a computational method needs to be able to account for.One example is so called spinpolarisation, which accounts for the fact that up and down spin electrons will respond differently to an excess of spin density, which in turn affects the spatial distribution of up and down spin electrons, meaning that it is no longer equal for each as is the case for a closed shell wavefunction.The wavefunction should also be constructed so that spin angular momentum is accounted for properly, i.e., the wavefunction should be an eigenfunction of the spin angular momentum operators describing the total spin (S) and the projection of spin along one direction (the operators Ŝ2 and Ŝz with eigenvalues S(S+1) and M s respectively).This ensures that the wavefunction describes a pure spin state, without being contaminated by terms from other spin states that describe additional unpaired electrons in the model, a problem referred to as spin-contamination.(Jensen, 2007) Highspin open-shell states such as X 3  g  can be treated qualitatively using a singleconfigurational model, while low-spin open-shell states (including a 1  g ) cannot.If one wishes to use a single configurational model then both the polarisation and contamination cannot be simultaneously accounted for.One can either use a restricted approach in which spin up and spin down electrons are maximally paired and share a common spatial component, or one can use an unrestricted approach in which each up and down spin electron have different spatial orbitals.In the restricted approach the single configuration is automatically an eigenfunction of Ŝz with an eigenvalue corresponding to half the excess number of spin up electrons (using the convention of the open-shell electrons having spin up), and the eigenvalue of Ŝ2 , corresponding to a pure value of S, given by the sum of spin up quantum numbers for all open-shell electrons.Since these open-shell electrons see the same spatial distribution of the paired closed shell spin up and spin down electrons there is no polarisation leading to non-physical effects such as no hyperfine coupling of excess spindensity to nuclear spin.On the other hand using an unrestricted approach allows for spinpolarization as the spin up and spin down orbitals can have different amplitudes in different regions of space.Such unrestricted wavefunctions mix in states of higher S (but common M s ).In general the higher this spin-contamination (i.e., how far above S(S+1) the  (Bally & Borden, 2007) Multi-configurational methods, as described in 2.1.3,are therefore essential for a balanced treatment of all the states in the oxygen molecule, and this includes when oxygen weakly interacts with a sensitizer, the topic covered in the next section.

Photosensitisation
In the discussion of the various aspects of the photochemistry that underpins the pathways involved in PDT, we have now arrived at the crucial stage where the reactive singlet oxygen species is generated through photosensitisation by the excited state of the irradiated photosensitiser molecule.The direct absorption of pure gaseous O 2 is very weak due to the nature of the electronic transitions, as discussed above, and in order to describe the generation of singlet oxygen through sensitisation, the nature of the perturbed O 2 molecule needs to be considered.(Badger, Wright et al., 1965) As mentioned previously, the perturbation is normally viewed as taking place through a collision complex between oxygen and the perturbing molecule, in our case the photosensitiser molecule, PS.This complex, which we will denote as PS-O 2 , and a selection of the lowest lying electronic states in the complex, is described in figure 11.
One of the first attempts to describe the mechanisms behind photosensitisation originated from observations involving heteroaromatic solvent molecules as perturbers of the O 2 molecule.The fact that oxygen interacts with heteroaromatic organic molecules was first inferred in the 1950s when it was shown that samples where oxygen had been dissolved gave rise to an absorption band at longer wavelengths than the aromatic compounds alone.
The first suggestion as to the nature of this interaction, and its precise mechanism, was presented in a seminal paper in the 1960s when Mulliken and Tsubomura suggested it to be due to absorption by a charge transfer complex (denoted 1,3 (PS + O 2 -) in Figure 11) originating from the interaction between oxygen and the aromatic perturbing molecule, where the oxygen has the role as electron acceptor.The interaction in these complexes were considered Fig. 11.Diagram illustrating some states of the PS-O2 complex.As new wavefunctions are needed to describe the complex, transitions that were forbidden in the isolated molecules can now become allowed.
to be very small, as the added electron in O 2  will vacate one of the already occupied π g orbitals, which have a very small overlap with the heteroaromatic n-donor orbitals.(Tsubomura, 1960) The mechanism of the enhancement of a transition in O 2 through a collision complex has been suggested to be due to a 'borrowing' of intensity from the spin allowed transition.The degree of enhancement in these cases will however depend on the nature as well as the orientation of the collider, as the degenerate anti-bonding π g,x and π g,y will be perturbed differently depending on the nature and the orientation of the valence MOs of the colliding sensitiser molecule.In the case of a photosensitiser the size of a porphyrin derivative, this is a crucial point in the modelling of the charge transfer event.
Several more possible mechanisms were suggested through the years, where the main theme was the enhanced intersystem crossing of the photosensitiser, where the states are mixing either directly or through charge transfer (CT) states of the complex.Kawaoka however later suggested that there might not be a need to populate the actual CT states of the PS-O2 complex in order to generate singlet oxygen, as its states were mixed to a varying degree.(Kawaoka, 1966) It was suggested that CT character will in fact have a negative effect on the production of O 2 (a 1 Δ g ).This holds true for intramolecular CT within the photosensitiser, which facilitates coupling of an excited state to a ground state, as well as for intermolecular CT.The latter is normally described as the degree of mixing of the CT states with the excited or ground state of the exciplex.(Kawaoka, 1966) If the PS-O2 complex has large CT character, radicals are easily formed, which has a negative effect on the singlet O 2 yield.Also, the CT state will facilitate coupling between excited PS-O2 states and the ground triplet state 3 ( 1 PS-3 O2).A further mechanism has also been suggested as it has been shown that O 2 ( 1 Δ g ) can be produced through a direct cooperative absorption by the PS-O2 ( 3 Σ g -) complex.(Pagani, Salice et al., 2010;Scurlock & Ogilby, 1988) This highlights the fact that many uncertainties still exist with regards to the understanding of the actual mechanism of the charge transfer pathway taking in the photosensitisation process.
It has become clear that O 2 quenches excited states of organic molecules of both multiplicities, with varying rates depending on the environment.The quenching of the S 1 state does not generally contribute to the generation of singlet oxygen, but leads to the formation of T 1 of the sensitiser, through O 2 enhanced ISC, or the generation of the sensitiser singlet ground state, O 2 enhanced IC, with regeneration of the O 2 ( 1 Σ g + ).(Schmidt & Schweitzer, 2003) With respect to PDT applications, quenching of the S 1 of the sensitiser can therefore be looked upon largely as a loss process.
In terms of modelling the interaction between the oxygen and the sensitizer one is faced with several problems.Firstly the weak nature of interaction complex, and the possibility of many geometrical conformers contributing in the condensed phase has to be considered.Secondly, and the open-shell nature of the collision complex requires the same techniques as isolated O 2 requires, i.e., multi-configurational methods.Given the often large size of the sensitiser this presents many challenges to theory.In fact it is fair to say that accurate multistate, multi-configurational investigations of this system have not been attempted to date.Modern methods (for example as discussed in section 2.2.3) can in principle be applied to the full interaction complex.These when combined with molecular dynamics simulations, and sampling, offer up the possibility that this final crucial step in PDT can be modelled as accurately as the other steps and thus hopefully lead to predictive use of modelling in the design of new PDT drugs.

Conclusions
We have introduced Photodynamic Therapy in terms of the photochemical pathways taking place within a target biological system.We presented an outline of the challenges associated with the modelling photosensitiser molecules, focusing on the issues related to electron correlation problem, in terms of methods using the electronic wavefunction as well as methods based on the electronic density of the system.Further we have highlighted the computational advances that have made it possible to accurately describe both linear and non-linear absorption processes in molecules the size of common photosensitisers.An accurate description of these processes can be crucial in the development of photosensitisers with absorption maximum within the optical window of tissue penetration.The modelling of ultrafast processes taking place from an excited state of a molecule is challenging and we have summarised the issues associated with the description of these processes in large molecules.Finally, we have discussed the difficulty in an accurate theoretical description of the energy transfer step involving molecular oxygen present in the target cell, mainly in terms of the challenges associated with the modelling of open-shell systems.
Despite several drawbacks, including poor light penetration of cancer tissue and poor targeting of the photosensitiser systems, PDT remains a very attractive goal, and research continues through many avenues.The search for the ideal photosensitiser, which will make PDT a realised anti-cancer treatment for general use is very much under way and the input that can be gained from computational studies cannot be ignored, as it has the potential to increase the understanding of the entire photochemical pathways involved.Computational techniques have now matured so that molecules the size of standard photosensitisers can be investigated routinely for their absorption properties.This is also true for non-linear absorption properties, which means that this approach to PDT can be developed and hence potentially improve the scope of the technique, and be an important tool both when rationalising the existing classes of photosensitizers and the design of future photosensitisers with desired molecular properties.
Even though significant challenges clearly still remain, the insight a theoretical input can give to the understanding of the fundamental photochemical mechanisms has potential to greatly accelerate the development of PDT as an everyday technique.

Fig. 4 .
Fig. 4. Summary of the TPA spectra.Note the change in scale between the two TPA crosssection axes.

Fig. 6 .
Fig. 6.Overview of the RASSCF strategy, where the active space is divided into three parts, RAS 1-3.

Fig. 7 .
Fig.7.Overview of the orbitals chosen in the description of the porphyrin chromophore, where the whole π-system have to be accounted for.RAS 2 includes the frontier molecular orbitals, FMOs, of the porphyrin macrocycle, as these have been shown to be vital for an accurate description of the absorption in molecules of this type.

Fig. 9 .
Fig. 9.The Porphyrin free base, with crucial structural sites indicated.All modern photosensitisers used in practice today are based on the porphyrin macrocyclic structure, and PDT with porphyrin derivatives have been under development since the 1960s (Fig 9).(Peng & Moan, 2003) In the body, the most recognisable porphyrin derivative that exists naturally is coordinated to Fe in the haem structure, present in the haemoglobin protein.The unchelated porphyrin structure is known as protoporphyrin IX (PpIX), which has been shown to be an effective photosensitiser.(Kennedy, Pottier et al., 1990) The strategy for PpIX treatment is to administer a synthetic version of δ-Aminolaevulic acid, the biosynthetic precursor to PpIX, which leads to a build-up of PpIX in the targeted area.In Europe, PpIX is approved for use in PDT of actinic keratosis and basal cell carcinoma.(Morton, McKenna et al., 2008) The most common photosensitiser in use today, however, is still the heamatoporphyrin derivative (HpD), which consists of a mixture of oligomeric porphyrin chains.In its purified form it is known as Photofrin®, which was the first PDT sensitiser to be approved for clinical use.(Elsaie, Choudhary et al., 2009) Photofrin or HpD is approved in Europe, the USA, Canada and Japan and is used mainly for cervical, bladder and gastric cancers.
expectation value  Ŝ2  U ) the worse the description of the state, and properties, including molecular geometries, are less accurate.Thus, one needs to go beyond a single configuration to describe these high-spin open-shell states more accurately, and such methods are the only starting point for the low-spin open-shell states.For a fuller description see the excellent review by Bally et al. which has a chemical discussion of all these aspects.