Biomolecular structure manipulation using tailored electromagnetic radiation: a proof of concept on a simplified model of the active site of bacterial DNA topoisomerase

We report on the viability of breaking selected bonds in biological systems using tailored electromagnetic radiation. We first demonstrate, by performing large-scale simulations, that pulsed electric fields cannot produce selective bond breaking. Then, we present a theoretical framework for describing selective energy concentration on particular bonds of biomolecules upon application of tailored electromagnetic radiation. The theory is based on the mapping of biomolecules to a set of coupled harmonic oscillators and on optimal control schemes to describe optimization of temporal shape, the phase and polarization of the external radiation. We have applied this theory to demonstrate the possibility of selective bond breaking in the active site of bacterial DNA topoisomerase. For this purpose, we have focused on a model that was built based on a case study. Results are given as a proof of concept.


I. Introduction
2][3][4] Similar methods, however, cannot be directly applied to biomolecules due to their complex structure.Moreover, the same biomolecule can contain a number of bonds like C-C or C-N, making it difficult to break them only on particular sites and not on others.
Particularly interesting biomolecules are enzymes like DNA topoisomerases, which solve the topological problems associated with DNA replication, transcription, recombination, and chromatin remodeling.Topoisomerases can be classified as type I and type II classes.This classification depends on whether they transiently break, respectively, one or two strands of DNA.Both type I and type II topoisomerases can be further subdivided into 5 subtypes: IA, IB, IC, IIA, and IIB. 5 One of the essential roles of topoisomerases is the strand-breakage reaction of DNA. 6Much interest has covered these enzymes, since they are currently important therapeutic targets involving numerous antimicrobial and anticancer treatments, and are also among the most effective targets. 7he building blocks of enzymes are amino acids, which exhibit large permanent dipole moments. 8When an external electric field is applied to a biomolecule, the dipoles tend to align parallel to the field.Since this cannot occur without substantial structural changes, an interesting interplay between electrostatic and conformational energies takes place. 9For this reason pulsed electric fields are used to modify the structures of proteins for industrial purposes.Depending on the intensity of the applied electric field, changes in the structure of the protein may occur which lead to malfunction.Therefore, pulsed electric fields (PEFs) acting on enzymes or proteins in general may produce structural transformation and affect their function.However, the field intensities needed for irreversible conformal changes might be too high and the applied fields might also destroy other structures.A possible way to avoid such collateral damage is to address and break particular bonds of proteins, leaving the rest of the system almost unchanged.
In this paper, we demonstrate that pulsed electric fields are not able to induce selective bond breaking and rather produce global conformational changes and damage.Therefore, a completely different approach is needed to address specific bonds.As mentioned above, many decades of intensive research on the field of ''selective light induced bond breaking'' have led to the possibility of breaking particular chemical bonds in small systems using electromagnetic radiation, but direct application to biomolecules is not possible.Different problems make selective bond breaking in biomolecules tremendously difficult.The most important one is the lack of specificity of simple energy inputs.Breaking C-C or C-N bonds in only one unit of a particular biomolecule not changing either the rest of the molecule or other molecules in the cell is an ambitious and challenging task, which cannot be achieved by merely applying light pulses with frequencies resonant to those of such bonds.The only promising approach is to ''prepare'' the electromagnetic pulses in such a way that they also contain information on the environment of the bonds.
Here, we present a theory for selective energy concentration and subsequent breaking of particular bonds of DNA topoisomerase.Our approach is based on electromagnetic-pulse shaping using optimal control (OC) schemes. 10The idea behind OC is the optimization of electromagnetic radiation (pulse shape, phase and polarization) in order to achieve a particular target state.For small molecules, optimal control of mode-selective femtochemistry has been described using quantum chemical molecular dynamics 11 We have applied our theory to simulate excitation of a model molecule consisting of 16 atoms, which represents the active site of topoisomerase.As a complementary result we performed an extensive analysis of the vibrational frequencies of the active site of DNA topoisomerase.Our results demonstrate the viability of selective bond breaking in biological systems using tailored electromagnetic radiation and constitute a proof of concept.

II. Modelling of topoisomerase
Three different models, each one including the catalytic tyrosine residue of topoisomerase III from Escherichia coli, were employed throughout the work: Model 1, which represented our base model for frequency determination (see Fig. 1); Model 2, a second model composed of 3 residues -the catalytically relevant tyrosine and the two immediate residues (see Fig. 1), and Model 3, obtained by adding the enzyme's environment and a single-stranded DNA (ssDNA) chain around Models 1 and 2 (i.e. the whole crystallographic enzyme).
The 3 models represent a growing sequence of the E. coli's topoisomerase III active site (Fig. 1); they were chosen to assess the influence of this growing complexity on relevant vibrational modes of the peptide structure (Model 1 and Model 2); to determine the eigenfrequencies and eigenvectors to be used as input for the optimization of electromagnetic radiation in order to achieve selective bond breaking (Model 1), and to assess the effects of external electric fields on protein conformation (Model 3).The models were built based on the crystallographic structure deposited in the protein data bank with code 1I7D. 12he phenylalanine residue (granting no enzymatic activity) was mutated back to a tyrosine (catalytically active).

III. Impossibility of selective bond breaking using constant electric fields
To demonstrate that pulsed electric fields are not able to break particular bonds avoiding other structural changes in biomolecules, we performed molecular dynamics (MD) simulations based on the GROMACS package. 13A total of 10 ns production dynamics were performed and additional 10 ns with an electric field dynamics were conducted.Details of the simulations are provided in the ESI † of this article.A detailed analysis was performed, regarding the effects of an external electric field on the protein structure of topoisomerase III.Three electric field settings were employed: 0.1, 0.5, and 1.0 V nm À1 .Notice that the magnitude of the external fields in our simulations should not be mistaken with the electric fields used at the industrial scale, 14 obtained by dividing the potential difference on a pair of electrodes by the distance between them.The electric field used at the industrial scale represents, in homogeneous systems, the net electric field that results from summing the external field arising from charging the electrodes with the counteracting field arising from the dielectric polarization of the medium.The actual effective electric field in the simulation medium can be determined by reducing the applied field by a factor equal to the relative permittivity of the medium. 15For example, given that the relative permittivity of TIP3P water is ca.90 16 applying an external field of 1.0 V nm À1 in our simulations produces an effective (net) field in the aqueous medium of ca.0.01 V nm À1 (100 kV cm À1 ), which is the same as the one used in, for instance, the food processing industry.
For electric field strengths of 0.5 and 1.0 V nm À1 relevant distortion of the enzyme's secondary structure was detected (see Fig. 2).We also observed significant changes in the rootmean-square deviation and radius of gyration, as well as in secondary structure elements for both field strengths, and in particularly for 1.0 V nm À1 (see ESI †).
These results provide us with the electric field ranges to be applied to all-atom simulations involving the full topoisomerase enzymes.What we observed is that the threshold value lies between 0.1 and 0.5 V nm À1 .What is also apparent is that with higher field strengths, the systems change their topology This journal is © the Owner Societies 2014 and assembly, faster.However, the field intensities needed for irreversible conformal changes might be too high and the applied fields might destroy other structures (collateral damage).
Additionally, to validate our protocol, a different protein from E. coli was selected: reduced thioredoxin (with the PDB code 1XOB). 17We applied the same protocol of the electric field dynamics to it.This was intended to evaluate the influence of the same electric fields on a totally different protein and to find out whether the same field strengths exerted the same influence on different proteins.We obtained for thioredoxin similar results to that for topoisomerase (see ESI †).From our simulations it is clear that application of pulsed electric fields near damage threshold can only destroy topoisomerase's active center and cause the unfolding of thioredoxin.This confirms the nonspecific character of a possible application of an external electric field to a certain medium and the possible effects of collateral damage.Therefore we conclude that pulsed electric fields cannot induce selective bond breaking in biomolecules.

IV. Our approach for selective bond breaking in biomolecules
When a molecule is irradiated with the help of an electromagnetic pulse, usually many vibrational eigenmodes are excited.This leads to oscillatory motion of all atoms.If the pulse is very intense, damage of the molecule through destruction of one or several chemical bonds occurs.The damage process is, in principle, uncontrollable if a single pulse is used.The goal of this paper is to show that it is possible to design a sequence of pulses able to induce selective bond breaking.The idea is sketched in Fig. 3.If the molecule is excited by a series of pulses, each one targeting a different eigenstate or groups of eigenstates, it might be possible to vary intensities, phases and eventually polarizations of the pulses in order to achieve breaking of a desired bond.
In order to investigate whether a shaped train of light pulses (as schematically shown in Fig. 3) is able to manipulate the bonds of a molecule, we need to start by achieving a proof of concept on a simplified model.In this paper we focus on the active site of topoisomerase and treat it as a gas-phase molecule.We then treat the molecule in the harmonic approximation in order to obtain its eigenstates.The interaction of electromagnetic radiation with the biomolecule is also modelled.Finally, we use an optimization scheme to obtain the proper parameters of the electromagnetic radiation able to concentrate such an amount of energy on particular bonds, so that they break.

IV.1 Mapping a biomolecule to a system of coupled harmonic oscillators
We consider a biomolecule consisting of N units (which could be atoms in an atomistic model or groups of atoms in the case of a coarse grained force field).The mass of each atom is M l (l = 1, N).In the harmonic approximation (see ESI †), the potential energy of the system is given as a function of the displacements u l of the atoms by where the indices l and n run from 1 to N while a and b refer to the directions in space (a, b = 1,3).The so-called force constants k lanb are the elements of the 3N Â 3N Hessian matrix K.Note that the real potential V eff ( u 1 ,. .., u N ) can be a complicated function of the displacements and is formally obtained after solving the electronic problem according to the Born-Oppenheimer approximation. 18It includes covalent, ionic and hydrogen bonds as well as van der Waals forces.In the limit of large relative displacement between two atoms the real (2)  where the vector -A(t) = (a 1 (t),a 2 (t),. ..,a 3NÀ6 (t)) contains the displacements written on the basis of the eigenvectors -G j of the dynamical matrix.O is the diagonal form of the dynamical matrix and contains the squared eigenfrequencies {x 1 2 ,. ..,x 3NÀ6 2 } of the molecule.

IV.2 Modeling the interaction of the molecule with electromagnetic pulses
Since an external electromagnetic field oscillating with frequency o will essentially excite the eigenstates of the molecule having the same frequency, we describe the displacements of the atom projected in the directions of the eigenvectors: a j ðtÞ Gj : In the presence of a sequence of n monochromatic pulses, the eigenvector displacements a j (t) will evolve (for all j) independently of each other according to the equations of motion where y(t À t i ) is the step (Heaviside) function, which ensures that the pulse is applied at t i and before this time no excitation with frequency o i takes place.E i is the product of the coupling strength and the amplitude of the ith pulse.This means, we consider the fact that the coupling strengths between the different modes and the external field are different.b j is the damping constant of the medium for mode j.Note that we will assume that the damping of all modes is the same and is determined by the viscosity of water, but the theory can also be applied if the damping for the different eigenmodes is different.n is the number of pulses applied.The initial conditions for the eigenvector displacements are a j (t j ) = 0 and : a j (t j ) = 0, since the molecule is at rest before excitation with electromagnetic radiation.

IV.3 Optimization of the electromagnetic radiation
The problem of pulse shaping consists of designing the best electromagnetic pulse in order to induce (favor) a particular reaction and at the same time to prevent other outcomes.This is a typical global search problem which relies on the maximization of a so called ''fitness functional'', which is in general a measure of the overlap between the actual molecular conformation at the control time and the target conformation.Different pulse shapes lead to different outcomes.The pulse shape maximizing the fitness functional is the optimal one.Now, we design a sequence of pulses in such a way that the frequencies o i of the electromagnetic waves are resonant with the vibrational eigenfrequencies o j .The optimization of the electromagnetic pulses will be based on the integration of eqn (4).We describe below an analytical and a numerical approach to optimize external radiation in order to achieve bond breaking.
Analytical optimization.In the resonant approximation, eqn (4) is reduced to With this approximation, which assumes that the molecule vibrates in its eigenstates and only resonant radiation can excite them efficiently, all equations of motion under the action of the electromagnetic radiation become fully decoupled and can be easily solved through direct analytical integration.Taking into account the initial conditions mentioned below (molecules at rest) the solution for t 4 t j ; is given by This solution holds when all vibrational eigenfrequencies are larger than the damping constant, which is valid for biomolecules.From eqn (6) we notice that the smaller the damping constant b j , the larger the desired final amplitude of a j (t) which is given by E j 2o j b j . Note also that the smaller the damping constant, the longer it takes for the system to reach the final amplitude, since the first two terms are only diminished Now we proceed to optimize the light pulse by using physical arguments.Our experience in optimal control theory 19,20 tells us that the best way of performing light control of atomic motion is to drive the system to the final target state as fast as possible.Therefore, we propose a fitness function of the form which becomes minimized if we choose the intervals between t final and the switch-on times t j to be the smallest possible time interval taking into account the system and the initial conditions.
It is important to recall that t final is the time at which the target displacements -A target are reached.In order to achieve concentration and eventually bond breaking, we have to choose -A target correspondingly.The phases of the pulses are determined by our choice of smallest possible t final we choose the switch-on time t j of each pulse in such a way that the desired amplitude a j (t final ) = a j,target is reached exactly after one oscillation period.This means that the last term of eqn (6), , must perform only one oscillation between t = t j (switch-on time for pulse j) and t = t final .This leads to the condition o j (t final À t j ) = 2p, or Eqn ( 8) reflects the fact that the pulse to drive the eigenmode j to the desired target amplitude a j , target must be switched This journal is © the Owner Societies 2014 on a time 2p/o j before the time t final at which the overall target displacements should be achieved.In order to simplify the expressions we initialize the time by setting t final = 0.Then, for each j it must hold a j (0) = a j,target .Replacing this condition into eqn (6) allows us to univocally determine the values of E j , which are given by Eqn ( 9) is valid for all j = 1, 3N À 6.In eqn (9), t j ¼ À 2p o j .
Eqn ( 8) and (9) determine the optimized (designed) parameters E j and t j for the electromagnetic pulses, so that the target displacements are reached at t = 0.
Thus, we were able to minimize analytically the fitness function of eqn (7) without using numerical approaches.Note that this is only possible in the resonant approximation for the coupling between the molecule and electromagnetic radiation, and within the harmonic approximation we used.Of course, if we go beyond the harmonic approximation and consider the full interaction of the molecule with electromagnetic radiation, the solution cannot be found any more analytically.However, the theoretical scheme presented in this work remains unchanged.
Numerical optimization.A major drawback of the fitness function F({E j }) of eqn (7) is that it only depends on the displacements at the final time t final and not on the whole time evolution of the bonds.This means that, even if the target displacements at the final time -A target are constructed in such a way that only one bond is significantly elongated while the others remain intact, F({E j }) cannot ''control'' whether the other bond elongate considerably on times smaller than t final .We correct now this problem and pay the price of a numerical optimization.
Note that we can integrate eqn (4) analytically without making the resonant approximation of eqn (5).Moreover, we can generalize eqn (4) to the case of the arbitrary phase of the electromagnetic radiation We can now construct a more complex and reliable fitness function, as follows.First, we define the relative timedependent elongation of the bond k between atoms j and l as l to its equilibrium position in the native structure.The index k labels the bonds.We relate d bond (k,t) to its value in the target by defining D k = d bond (k,t)/d bond (k,target).D k = 1 in the ideal case of perfect match of bond k with its length in the target structure.Now, the fitness function must be such that a particular bond becomes dramatically stretched at the target time, while all others remain only slightly elongated at all times.Only in this way selective energy concentration on one bond can be guaranteed.This requirement is much more difficult to fulfil than that imposed by the previous fitness function of eqn (7).Therefore, we define the new fitness function J as y(x) refers to the step function, and the factor e tells us how much the rest of the bonds may be elongated with respect to a particular target situation.Note that in this way all bonds different from the target bond will have at most an elongation equal to 100e% of the target bond.We choose e = 0.55.
Bond breaking vs. energy concentration.It is important to point out that in the harmonic approximation only light induced concentration of energy can be achieved.Bond breaking cannot be described in the harmonic approximation.However, it is possible to estimate if bond breaking can occur.One can construct the fitness function by requiring that the light pulse to excite a particular combination of eigenstates, so that only the red and dark blue atoms (see Fig. 4) move with a kinetic energy larger than the dissociation energy, i.e. sufficient to break the bond, while the other atoms only perform small amplitude oscillations with low kinetic energy.

IV.4 Frequency analysis protocol
In order to find out the optimized electromagnetic radiation through application of the theory presented in the previous subsections, we had first to determine the vibrational eigenfrequencies and eigenvectors of the biomolecule.For this purpose we used the harmonic approximation described in the ESI.† Our study was directed to finding out the vibrational frequencies of a system with few amino acids.This is important since the Fig. 4 Schematic description of the interaction of an electromagnetic pulse with a system of coupled harmonic oscillators (which represents part of a biomolecule).The optimally shaped pulse excites the system in such a way that it induces such large amplitude oscillations on atoms belonging to a particular bond (marked by a red dashed line), that it breaks, while all other atoms only perform small oscillations.Bond breaking means that the kinetic energy of the red and dark blue atoms involved in the bond is larger than the dissociation limit, as is illustrated in the sketch of the potential corresponding to the mentioned bond.vibrational frequencies are highly influenced by the environment, namely for proteins due to the strong polar topology that they display.In addition, their intermolecular interactions dramatically influence the vibrational dynamics and structural properties of vibrational frequencies. 21However, one cannot easily determine the vibrational frequencies of an entire protein, since it requires a tremendous computational cost.Due to this aspect, we have started from the simplified Models 1 and 2 described above.The frequencies of selected vibrational stretching modes were determined for each of the above models, in which we focus our analysis towards the backbone atoms, including the C a and C b carbon atoms, as well as the atoms of the trans peptide group.They represent important bonds on the peptide structure.Both gas-phase and a water solvation model were used for frequency calculations.For the solvation model we have used the Self Consistent Reaction Field (SCRF) theory, implemented through a Polarized Continuum Model (namely the Conductor-like Polarizable Continuum Model, i.e.CPCM). 22The B3LYP functional, 23 within the Density Functional Theory (DFT), coupled to a 6-31G(d) basis set 24 was employed since good results are reported for this combination when considering frequency calculations, and also due to the small computational cost. 25he frequencies were multiplied by a scaling factor of 0.9614 as is customary in frequency calculations with this functional and basis set. 26This is necessary since ab initio harmonic vibrational frequencies are typically larger than experimental fundamentals.This scaling can also be applied due to the uniformity of this variation between theory and experiment.
The same chemical bond vibrational frequencies were compared in both models, in order to assess the biomolecule-specificity of a radiation-induced bond breaking event in the protein backbone.In addition, the influence of an external electric field (1.0 V nm À1 ) was accessed for the variation in vibrational frequencies, in Model 1 and Model 2. All calculations were performed using Gaussian 09. 27

V. Results and discussion
In this section we demonstrate that it is possible to induce energy concentration on a particular predetermined bond of Model 1 for the active site of topoisomerase by applying an optimized train of light pulses.We first show results corresponding to the frequency analysis protocol for Model 1 and Model 2, obtained from the computation of the mass weighted Hessian matrix and its diagonalization.The so obtained eigenvectors and eigenvalues for Model 1 are then used as input for the pulse tailoring.

V.1 Frequency analyses
We have calculated the vibrational frequencies in Models 1 and 2 to later study selective bond breaking in biomolecules.It is important to study both models because vibrational frequencies should be sensitive to the environment.It has been described that vibrational modes of proteins are affected by the protein's structure and sequence, the exposure of the surface residues to the solvent (mediated through hydrogen bonding), and by bulk effects. 21Other factors can influence the vibrational modes of proteins as for instance the protonation states of amino acid residues or the interaction with ions.This sensitivity to the chemical composition and the architecture of molecules renders infrared spectroscopy as one of the classical methods for investigating protein structure, protein-mediated reactions, and protein folding.
In the ESI † we present a detailed summary of all the results; here we present only the relevant conclusions to avoid a lengthy, repetitive description of the same kind of results.
We have focused our analysis in the backbone atoms represented in our models, specifically the trans peptide group and the C a and C b atoms.The results are discussed relatively to the gas-phase frequencies; however the water dielectric determined frequencies are also compared, since a relevant deviation between the two was observed for specific vibrational modes.The trans peptide group vibrations include amide A, and amide I, II, and III vibrations.
The relevant vibrations will be further discussed, compared between both models, and also with what is expected experimentally.Considering Model 1, the results show a relevant difference between the gas-phase and the water dielectric calculations, specifically for those vibrations within the range 1709-1741 cm À1 and for low frequency modes.The frequencies between 1709-1741 cm À1 can be related to amide I vibrations of the simplified peptidic model.This difference can reach ca.30 cm À1 .For Model 2 (zwitterion) the results are similar.We also observe high deviation within the 1658-1708 cm À1 range, which also corresponds to the amide I vibration.This difference for the amide I vibration has been previously characterized in the literature. 24This supports the influence of the environment on the vibrational frequencies of amino acids.Obviously, more thorough characterization by representation of the protein environment would be necessary to attain a more thorough description of the influence of the specific environment on the studied frequencies, but this is out of the scope of this paper.
When comparing the same vibrations on both Model 1 and Model 2, we observe relevant shifts in the analyzed frequencies.Some relevant observations of these shifts, focusing on amide A, I, II, and III vibrations, are depicted in Table 1.Additionally, vibrations involving the bonds of interest are depicted also in Table 1 (C a -C b , C a -C bb , and C a -N bb stretching modes).
Amide A vibrations correspond to the NH stretching vibrations, normally found between 3270-3310 cm À1 .
For Model 1 we see this stretching event at 3466 cm À1 , whereas for Model 2 this is found at 3437 cm À1 .Amide I vibrations (near 1650 cm À1 ), in general, have a high CQO stretching component and minor contributions from the out-of-phase CN stretching, the CCN deformation and the NH in-plane bending.In Model 1 this mode is found at 1709 cm À1 and in Model 2 at 1700 cm À1 .A second amide I contribution is found at 1741 cm À1 and 1625 cm À1 , for Model 1 and Model 2, respectively.The discrepancy between both models must be related to the to the change of the immediate environment.In Model 2, this second amide atom is attached to two carbon atoms (proline residue).The amide II stretching (normally found at ca. 1550 cm À1 ) is the This journal is © the Owner Societies 2014 out-of-phase combination of the NH in plane bending and the CN stretching vibration with smaller contributions from the CO in plane bending and the CC and NC stretching vibrations.For Model 1 this is found at 1481 cm À1 and for Model 2 at 1501 cm À1 .N-deuteration converts the amide II mode largely to a CN stretching vibration at 1460-1490cm À1 (amide II 0 ).For Model 2, this CN vibration is in the range 1416-1426 cm À1 , and it is probably related to the N atoms belonging to a proline residue.
Other interesting modes are related to the C a -C b stretching: at 1386 and 1074 cm À1 for Model 1, and at 1349 and 1007 cm À1 for Model 2. The presence of C a -N bb stretching, at 1186 and 1040 cm À1 (Model 1 and Model 2, respectively), is also relevant.Finally, the C a -C bb vibration is also found for Model 1 at 1355 cm À1 , and for Model 2 at 1416 cm À1 .
Our results confirm that the vibrational frequency of a given bond changes upon changing the immediate environment (here represented by the increase of the molecular model).This corroborates the fact that vibrational frequencies are highly influenced by the environment due to the strong anisotropic polar topology that proteins display.Also intermolecular interactions influence the vibrational dynamics of the chemical bonds.

V.2 Proof of concept: selective bond breaking using tailored electromagnetic radiation
We show in the following that it is possible to optimize the shape and the phase of electromagnetic pulses in order to induce selective breaking of one bond in Model 1 of topoisomerase, while preserving all other bonds.
The goal of this work was to tailor pulse trains able to break three bonds of the molecule, as shown in Fig. 5.
As an example we show in this section the analytical optimization of the pulse train 3 which must be able to break only the bond between C a and C b (Fig. 5) in Model 1 of topoisomerase.Analytical results for the optimization of pulses 1 and 2 are included in the ESI.† Moreover, we show here also results of the more accurate numerical optimization for pulses 1 and 3. Pulse 1 is designed to be able to break the bond between N bb and C a (Fig. 5).
Results of the analytical optimization.We first mapped Model 1 of the active site of topoisomerase to a system of coupled harmonic oscillators by applying the harmonic approximation, as described in Section IV.From the computations on the Hessian matrix described therein we extracted the vibrational eigenfrequencies o j ( j = 1, 3N À 6) and the eigenvectors -G j ( j = 1, 3N À 6).With this information, the unitary transformation matrix C was constructed (see ESI †).Also the matrices M, M 1 2 and M À 1 2 were computed from the knowledge of the masses of the atoms composing the Model 1 segment.
As a first output, the switch-on times t j of the 42 pulses were determined using eqn (S5) of the ESI † and taking into account of that we set t final = 0. Since we apply resonant electromagnetic pulses, in principle 3N À 6 pulses are needed.Since Model 1 of topoisomerase consists of 16 atoms, the designed trains contain, at most, 42 pulses.
In order to obtain the values E j of the electromagnetic pulses using eqn (9), we constructed first the target-displacement vectors -U target in real for the three different final states (energy concentration in each of the three mentioned bonds).Finally, we project them into the eigenvectors to obtain -A target using the relationship In Fig. 6 we show the form of the target-displacement vectors to which pulse train 3 was designed to drive (see Fig. 5).Notice that we constructed the target vector in such a way that the bond to be broken acquires the largest possible stress while the other bonds perform at most small oscillations at the final time t = 0.
Once the vectors -A target are obtained from the real-space target displacements constructed in Fig. 6 and 7, they are inserted in eqn (9) to determine the amplitude of the pulses in the trains.This is the solution of the problem.To demonstrate that this solution works we put those amplitudes and the previously determined switch-on times in eqn (6) and obtained the time behavior of the eigenvector-displacements. Then we converted them to real displacements and calculated the time evolution of all 15 bonds of the molecule.The performed determination of the tailored pulse trains has and the masses as input.As output we obtain t j ( j = 1, 3N À 6), E j ( j = 1, 3N À 6) a j (t), and the bond lengths vs. time.Fig. 5 For Model 1 of topoisomerase we design three pulse trains able to selectively break the three marked bonds.Labels 1, 5, 7 and 10 correspond to carbon atoms (colored grey); 2 and 11 to oxygen atoms (colored red), and 3 and 12 to nitrogen atoms (colored blue).Pulse train 1 will be designed to concentrate energy in the bond between atoms 3 (nitrogen) and 5 (carbon), pulse train 2 will break the bond between atoms 10 (carbon) and 5 (carbon) and finally pulse train 3 will induce concentration of energy and bond breaking between atoms 7 (carbon) and 5 (carbon).
Fig. 7 shows the effect of the designed pulse train 3 on Model 1 for topoisomerase.Note that the bond between C a and C b performs strong oscillations, while other bonds remain with fluctions under 25%.It has to be mentioned that a bond fluctuation of 20% is very likely to be enough to break the bond.However, we can scale down all bond fluctuations by requiring the maximum bond length at the target time to be 50%, for example.Taking into account that the optimization has been performed fully analytically and using a very simple fitness function (eqn (7)), this is a very good starting point and shows that the concept works.We show below that the improved fitness function together with its numerical treatment successfully avoids large fluctuations in the bonds that are not selected to be broken.
Results of the numerical optimization.We solved the equations of motion of eqn (10) analytically and then proceeded to perform numerical optimization of the fitness function J({E j ,j j }) eqn (11) using the SLSQP optimizer, 28 which is a sequential least squares programming algorithm based on the Han-Powell quasi-Newton method.We used a friction parameter b = 47 cm À1 obtained from the Literature. 29In order to rule out a bias in the optimization algorithm we have also numerically optimized the fitness function J({E j ,j j }) using the L-BFGS-B algorithm, 30 which basically gave the same results.Results of the relative bond elongations vs. time upon optimized pulse trains 1 and 3 as well as the shape, spectral density and phase of the optimized pulse trains are shown in Fig. 8-11.
Fig. 8 shows the behaviour of the bonds in Model 1 of topoisomerase as a function of time upon application of the pulse train 3, designed to break the bond between C a and C b (Fig. 5).Fig. 8 has to be compared to Fig. 7, since both show the response of Model 1 to differently optimize pulse trains designed to break the same bond.It is clear that the fitness function J({E j ,j j }) works much better than F({E j }), since it forces all not selected bonds to perform small oscillations during the whole time.Another important difference between the analytical and the numerical optimizations is given by the time-scale.Whereas the pulse trains designed analytically have a total width of several picoseconds, the numerical optimization is able to produce selective bond breaking on the femtosecond scale.However, the numerical optimization is not restricted to ultrashort time scales.The width of the pulse train can be varied in eqn (11).
Fig. 9 shows the numerically optimized pulse train 3. The pulse envelope has a smooth shape, while the spectral density indicates that frequencies resonant with certain vibrational modes of Model 1 have a large contribution.
Finally, in Fig. 10 and 11 we show results for the numerical optimization of pulse train 1 designed to break the bond between N bb and C a (Fig. 5).Again, all other bonds only perform

Fig. 1
Fig. 1 Representation of the three models used in this study.Hydrogen atoms were omitted for clarity.Model 1 is rendered in pink and Model 2 in pink and magenta.Model 3 represents the full protein, represented in silver and green and with the ssDNA colored in blue.

-
u N ) goes to zero.This indicates the dissociation limit.In the harmonic approximation [eqn (1)], however, V harm eff ( u 1 ,. .., u N ) increases parabolically for increasing relative displacements.If the biomolecule is in its native state (global minimum) or in a local minimum of V eff the displacements of the atoms (units) acquire their equilibrium values equal to zero.From eqn (1) one can easily derive the equations of motion for the displacements (see ESI †).By projecting the displacements into the directions of the eigenmodes of the molecule, which are determined from the diagonalization of the mass weighted Hessian (dynamical) matrix (see ESI †), one obtains € ÃðtÞ ¼ ÀO ÃðtÞ;

Fig. 2
Fig. 2 Molecular representations and results for each of the MD simulations conducted -at the end of 10 ns simulation in each situation, with the exception of the 1.0 V nm À1 situation, in which the structure was chosen until the periodic boundary conditions were kept (ca.1.4 ns).Only the protein structure (green and silver surface) and ssDNA chain (orange) are depicted.Water molecules are omitted for clarity.

Fig. 3
Fig. 3 Scheme of selective bond breaking through a train of electromagnetic pulses.Different monochromatic pulses of frequencies o i are applied at different times t i , which lead to bond breaking at the target time T target .
where -R j (t) refers to the time-dependent position of atom l and -R 0

Fig. 6
Fig. 6 Representation of the vector U ~target to be reached by Model 1 after excitation through the tailored pulse train 3. Arrows indicate the direction of the displacements.Displacements of atoms 5 and 7 are set antiparallel to each other and large in magnitude.The other displacements are chosen in such a way that the molecule can effectively be broken at the desired position.Thus, only the bond between atoms 5 and 7 is strongly stretched, while the other bonds do not suffer considerable stress.

Fig. 7
Fig. 7 Effect of the tailored pulse train 3 on the time dependence of the 15 bonds of Model 1 of topoisomerase.The black vertical lines represent the times at which the different pulses of the designed train are switched on.The red line shows the behavior of the bond length between atoms 5 and 7.Note that this bond is elongated by 100% at the final time, whereas all other bonds (green lines) are at most compressed/elongated by about 25%.

Fig. 8
Fig. 8 Time evolution of bond lengths in Model 1 of topoisomerase upon excitation with the numerically optimized pulse train 3. The blue line shows the bond length between atoms 5 and 7. Note that this bond is elongated by almost 30% at the final time, whereas all other bonds are at most elongated by less than 10%.

Fig. 9
Fig. 9 Numerically obtained optimal pulse train 3. (A) Shape of the pulse envelope in units of force over mass.(B) Spectral decomposition of the pulse train in terms of the vibrational frequencies of Model 1.The relative strength is calculated as E j /bo j .The phase is also indicated by the color of each peak according to the scale on the right.

Fig. 10
Fig. 10 Analogous to Fig. 8 for the response of bond lengths to the numerically optimized pulse train 1.The violet line shows the bond length between atoms 3 and 5.This bond is elongated by almost 25% at the final time, whereas all other bonds are at most elongated by 12%.

Fig. 11
Fig. 11 Numerically obtained optimal pulse train 1. (A) Shape of the pulse envelope.(B) Spectral decomposition of the pulse train in terms of the vibrational frequencies of Model 1.The phase is also indicated.

Table 1
Selected frequencies for relevant backbone vibrational modes and their respective differences in Model 1 and Model 2 a These modes are often coupled to other vibrational modes.n -bond stretching.