Fatima Zohra Benidris
email: fatima.benidriss@gmail.com

Badr Benmammar
email: badr.benmammar@gmail.com

Tarik Fethi

Bendimerad
email: ftbendimerad@gmail.com

Comparative studies of artificial intelligence techniques in the context of cognitive radio

Keywords: component: cognitive radio, cognition cycle, artificial intelligence

Cognitive radio (CR) is considered as the key enabling technology of Next Generation Wireless Systems (NGWS). In this context, CR enables users to dynamically access and to fairly share the spectrum with other users. Our paper describes some artificial intelligence techniques like artificial neural network, metahueristic algorithms and hidden Markov model, these techniques are proposed to provide the cognition capability in a cognitive engine.

INTRODUCTION

Wireless spectrum is a costly resource which is licensed by governmental agencies to operators for long periods of use. However, a large portion of the assigned spectrum is used sporadically leading to under-utilization of valuable frequency resources. To address this critical problem, Federal Communications Commission (FCC) has recently approved the use of unlicensed devices in licensed bands. Consequently, dynamic spectrum access techniques are proposed to solve these current spectrum inefficiency problems. This new area of research foresees the development of cognitive radio (CR) networks to further improve spectrum efficiency [START_REF] Wu | [END_REF]. The basic idea of CR networks is that the unlicensed devices (also called cognitive radio users) need to vacate the band once the licensed devices (also known as primary users) are detected. CR networks, however, impose a great challenge due to the high fluctuation in the available spectrum as well as diverse quality-of-service (QoS) requirements [START_REF] Wu | [END_REF]. In [Haykin, 2005], CR is an intelligent wireless communication system that is aware of its environment, can also learn from experience and can make changes to certain operating parameters (e.g. transmit-power, carrier-frequency, and modulation strategy) to adapt to the incoming RF stimuli in realtime. The main objective of a CR is to have a highly reliable communications with efficient utilization radio spectrum in order to satisfy the user needs [Baldo, 2008]. From the above definition, there are six key words stand out in this definition: aware, intelligent, learn, adapt, reliable and efficient. This paper will focus only on the learning ability of a CR. Cognitive radios will be employed according to a cognition cycle that was originally described by [Mitola, 2000] as the fundamental activities in order to interact to the environment. Figure 1 shows the activities that a CR should perform: observation, orientation to determine its importance, creates alternative plan, make decision and then implement the actions. Finally, is the learning activity that uses observation and the outcome of the decisions to improve the radio operation [Baldo, 2008]. These knowledge gathering will be exploited in the future orientation activities to produce a more effective decision.

Figure 1 : Mitola's simplified cognition cycle [START_REF] Ruslan | [END_REF].

Stated by [Rondeau, 2006a], CRs are all about optimization to suit the applications and user's needs.

It is common to practice the cognition cycle in solving optimization problem. Finding effective solution for an optimization problem in implementing CR is not an easy task. Cognitive engine which is the core of a CR will performs the modeling, learning and optimization process. It will reconfigure the communication system depending on which the radio operates [Rondeau, 2006b]. The engine must be able to interpret the user needs to the radio action. This continuous activity must involve either through user intervention or assuming that the user is expected to understand its needs from the radio. The CR engine should consider all the factors that it need to learn and make the best decision.

Correct learning data will produce correct decisions. The remainder of the paper is organized as follows Section two describes the cognitive radio architecture and discusses the reasoning and learning engines. Section three reviews AI techniques proposed for use in a CE and presents examples of their applications. Section four we take a look at the comparison of different AI techniques. Conclusions are drawn in Section 5.

COGNITIVE RADIO ARCHITECTURE

A software radio (SR) can be defined as a radio implemented with generic hardware that can be programmed to transmit and receive a variety of waveforms. Cognitive radio is often thought of as an extension to software radio, and here we treat it as such. A cognitive radio extends a software radio by adding an independent cognitive engine, composed of a knowledge-base, reasoning engine, and a learning engine, to drive software modifications. A well-defined API dictates communication between the cognitive engine and the SR. Figure 2 illustrates this architecture and the interaction between various components.

Figure 2 :

The architecture of the CR and the interaction between different components [Clancy, 2007].

Types of CR

According to [START_REF] Ruslan | [END_REF] there are three levels of CR systems which are basic CR, reasoning CR and learning CR.

A basic CR: is a radio system that senses and adapts itself to the environment but does not necessarily include any reasoning and learning techniques. It will measure all available channels and select one channel with the minimum noise and interference level. At any given time, the cognitive engine generates conclusions based on information defined in the knowledge-base, the radio's long-term memory. These conclusions are the result of extrapolations of this information based on reasoning or learning.

A reasoning CR: is an improvement of the basic CR with no interference that can automatically determines legal limits on secondary spectrum operation. The reasoning engine is often referred to in AI literature as an expert system. The learning CR: will always updates the decisionmaking based on the previous action taken. The learning engine is responsible for manipulating the knowledge-base from experience. As lessons are learned, the learning engine stores them in the knowledge-base for future reference by the reasoning engine. Depending on the application, the learning engine may only be run to train a newly initialized radio, or it could be run periodically as the radio operates. Standard AI techniques may be useful to use for the learning process in CR.

ARTIFICIAL INTELLIGENCE TECHNIQUES FOR COGNITIVE RADIO

This section presents some AI techniques that have been proposed throughout the literature as possible candidates for CR. They are presented in the order of historical development.

ANN (Artificial neural network)

The first artificial neural was presented by the neurophysiologist W. McCulloch and the logician W.

Pits in 1943 for the study of the human brain. The idea of artificial neural network (ANN) was then applied to computational models. Modeled on a nerve plexus, an ANN is nothing more than a set of non-linear functions with adjustable parameters to give a desired output [Haykin, 1999]. Different types of ANNs are separated by their network configurations and training methods, allowing for a multitude of applications. However, they are all comprised of neurons interconnected to form a network. Each artificial neuron usually produces a single output value by accumulating inputs from other neurons. While there are many types of ANNs available in the literature, only those most common and applicable to CR are presented here.

• Multi-layer linear perceptron networks: are comprised of layers of neurons, each being a linear combination of the previous layer's outputs.

• nonlinear perceptron networks although multilayer: can provide highly flexible and dynamic results, their network configuration must often reflect the data that they represent • Radial basis function networks: have a built-in distance criterion with respect to a center (a radial nonlinear function) in its hidden layer.

Figure 3 :

Typical example of ANN [START_REF] Ruslan | [END_REF].

Application of ANN to CRs: Because of their ability to dynamically adapt and be trained at any time, ANNs are able to "learn" patterns, features, and attributes of the system they describe. The term "learn" refers to the fact that the neurons are stored in computer memory, the outputs of which can systematically be adjusted to yield a new result for a new situation and remember the results. The attributes can be highly nonlinear, complex, and numerous, yet ANNs can be constructed by only a few examples, thus reducing the complexity of the solution.

For this reason, they have long been used to describe functions, processes, or classes that are otherwise difficult to analytically formulate. Therefore, ANNs can be used not only to classify or recognize received stimuli but to assist in the solution adaptation process as well.

The ANN has been adopted in spectrum sensing for CR [START_REF] Fehske | [END_REF]. The ANN has also been used for radio parameter adaptation in CR [Reed, 2005].

Baldo and Zorzi propose to use the ANN to characterize the real-time achievable communication performance in CR. In [Amraoui, 2011], the authors describes a new approach using CR to improve wireless link reliability for a cognitive radio mobile terminal; it is based on supervised learning. In addition, the ANN has been used for pattern classification in a pattern-based transmission for CR [START_REF] Ustundag | [END_REF], [START_REF] Orcay | [END_REF].

Metaheuristic Algorithms

Explicit relations between the parameters of a CR and the desired performance metrics are usually not available. Therefore, search algorithms based on mathematical relations cannot be applied to find the optimal parameters with respect to the performance metrics. Instead, metaheuristic algorithms [START_REF] Blum | [END_REF] can be applied to computationally hard problems to search through the solution space while learning and establishing the requisite relationships. Although the term "metaheuristic" was probably first mentioned in 1986 [Robbins, 1951], it can be traced back to earlier work on stochastic optimization methods in the 1950s [Neel, 2007].

Application of Metaheuristic Algorithms to CR:

The metaheuristic techniques presented here can not only be used for reasoning or finding the optimal solution with objective/utility function but can also be used for learning with the aid of training examples when the relationship between parameters and a desired performance measure is not well understood.

The objective of learning is to identify a hypothesis or a rule set from the search space that maximizes the fit of the training examples to the target concept or, in other words, to identify a hypothesis set or a set of rules that is consistent with the training examples.

Although the characteristics of each search algorithm are different, as can be seen in Table I, a common challenge in the application of metaheuristic techniques is the formulation of extensive examples for target scenarios. Among the various metaheuristic algorithms, the GA has been widely adopted to solve multi objective optimization problem and dynamically configure the CR in response to the changing wireless environment [Rondeau, 2004a].

HMM

The hidden Markov model (HMM) was first introduced in the late 1960s. It is a convenient and mathematically tractable statistical model to describe and analyze the dynamic behavior of a complex random phenomenon [Rabiner, 1989] that can be modeled as a Markov process with observable and unobservable states. The HMM generates sequences of observation symbols by making transitions from state to state, one symbol per transition. However, the states are hidden, and only the output is observable.

In general, a real-world process can be expressed as a random process producing a sequence of observation symbols or patterns with hidden parameters generating the observables. The symbols or patterns may be discrete or continuous depending on the specific processes. An HMM can be built for a specific system to explain and characterize the occurrence of the observed symbols or patterns. This model can then be used to identify the sequences of observations with the same pattern by choosing the model that would most likely produce the observed sequences. Therefore, an HMM can be used as an observation process of the CE to recognize or classify received stimuli and can achieve awareness. In addition, since it can reproduce the training sequences, it can be used for prediction. Furthermore, learning can be accomplished by creating new models. HMMs have been applied to CR research. Rondeau et al. propose to model the wireless channel online using an HMM for CR [Rondeau, 2004b]. The HMM is trained using the GA with data from a broadband channel sounder in a line-of-sight additive white Gaussian noise channel.

HMMs have also been used for spectrum sensing in CR [START_REF] Kim | [END_REF]. In addition, HMMs have been used for spectrum occupancy prediction [START_REF] Akbar | [END_REF].

COMPARISON BETWEEN THE DIFFERENT

TECHNIQUES

The table below shows the AI techniques that have been described in section three while showing the strengths and limitations and options of each technique, so the choice of one or more techniques at a time (the combination) is made according to user needs and the efficient use of spectrum.

Algorithm

IV. CONCLUSION

AI techniques lie at the heart of CR, and understanding the tradeoffs in the selection and design of AI processes is critical to a successful CR design. This paper has reviewed some AI techniques-ANNs, metahueristic algorithms, HMMs-that have been proposed to provide the cognition capability in a CE. While we have seen that AI techniques have been applied to numerous CR applications [Baldo, 2008]- [Glover, 1986], many implementations remain rudimentary, perhaps due to the interdisciplinary nature of the field and perhaps because products are just beginning to appear. We have seen that the appropriateness of AI techniques varied by application and implementation.

The decision in choosing one or some AI techniques over other techniques in CE design needs to be made based on the application requirement, considering the tradeoffs among response time, processing complexity, training sample availability, robustness, etc. In addition, the learning capability of the AI technique needs to be considered and exploited in designing a CE as learning is critical to the performance of autonomously deployed CRs.

Table 1 .

 1 Comparison between the different techniques

		Strengths	Limitations	Options
	Artificial neural network	Ability to describe a	Training may be slow	Can use other
	(ANN)	multitude of functions ;	depending on network	learning techniques in
		Conceptually easy to	size ;	the training phase
		scalable;	Possible over	(i.e.,GA) ;
		Excellent for	training ;	Can be combined
		classification ;	No theory to link	with RBS (rule-based
		Can identify new patterns ;	application with	system);
			required network ;
	Metaheuristic Algorithms Excellent for parameter	Formulation of rule	Can be used in
		optimization and learning	space is difficult when	conjunction with
		involving relationship	learning or	RBS ;
		between parameter values ;	optimization is not	Learning can also be
		Can use other learning	restricted to parameter	used in the search
		techniques in the training	values ;	process ;
		phase (GA).
	Hidden Markov model	Can model complicated	Requires good training	Based on previous
	(HMM)	statistical processes ;	sequence ;	knowledge, RBS et
		Good for classification ;	Computationally	CBS can help HMM
		Easily scalable ;	complex ;	determine the
		Can predict based on		observation duration
		experiences ;		for a specific
				application and
				overcome issues with
				new situations ;


B Tomas Johansson
email: b.t.johansson@bham.ac.uk

Daniel Lesnic

Thomas Reeve

A method of fundamental solutions for two-dimensional heat conduction

Keywords: Heat conduction, Method of fundamental solutions, Regularization, Two-dimensional Domains, Collocation heat conduction, method of fundamental solutions AMS Subject Classification: 35K05, 35A35, 65N35

 for one-dimensional heat conduction with the sources placed outside the space domain of interest, to the two-dimensional setting. Theoretical properties of the method, as well as numerical investigations, are included, showing that accurate results can be obtained efficiently with small computational cost.

Introduction

In many engineering applications finding the solution to various heat conduction problems is of fundamental importance. Examples include, heat exchangers, mathematical finance, in particular after transforming the Black-Scholes equation into the heat equation, and various chemical and biological systems, including diffusion and transportation problems. Thus, due to its importance, many different numerical techniques have been developed for calculating heat flow. The method of fundamental solutions (MFS) is a powerful numerical technique that has been used to obtain highly accurate numerical approximations of solutions to linear partial differential equations (PDEs) with small computational effort, see the reviews [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Golberg | The method of fundamental solutions for potential, Helmholtz and diffusion problems[END_REF]. However, this technique has mainly been applied to stationary heat flow governed by elliptic partial differential equations [START_REF] Balakrishnan | The method of fundamental solutions for linear diffusionreaction equations[END_REF][START_REF] Bogomolny | Fundamental solutions method for elliptic boundary value problems[END_REF]. Recently, in [START_REF] Johansson | A method of fundamental solutions for transient heat conduction[END_REF], an MFS for the time-dependent linear heat equation in one spatial dimension was proposed and investigated. This method was extended to free surface Stefan problems in [START_REF] Chantasiriwan | The method of fundamental solutions for free surface Stefan problems[END_REF] and to heat conduction in one-dimensional layered materials in [START_REF]A method of fundamental solutions for transient heat conduction in layered materials[END_REF]. Encouraged by these results, in this paper we extend the approach considered in [START_REF] Johansson | A method of fundamental solutions for transient heat conduction[END_REF] to heat conduction in two-dimensional bodies. We note that other formulations of the MFS for the parabolic heat equation were given in [START_REF] Chantasiriwan | Methods of fundamental solutions for time-dependent heat conduction problems[END_REF][START_REF] Golberg | The method of fundamental solutions for time-dependent problems[END_REF][START_REF] Kolodziej | Transient heat conduction by boundary collocation methods and FEM -A comparison study[END_REF][START_REF] Valtchev | A time-marching MFS scheme for heat conduction problems[END_REF][START_REF] Young | Time-dependent fundamental solutions for homogeneous diffusion problems[END_REF].

We begin the work in Section 2 by introducing some notation and function spaces, and formulate and review results for the linear heat conduction equation. In Section 3, we prove some theoretical results that we will use in our MFS formulation, in particular, that linear combinations of fundamental solutions are dense in the space of square integrable functions defined on the boundary of the solution domain. In Section 4, we describe the MFS implementation for various (bounded) solution domains. In Section 5, we present numerical results for both some direct problems and also for one inverse problem. These results show that accurate approximations can be obtained efficiently with small computational effort.

Preliminaries and notation

We let x = (x 1 , x 2 ) and y = (y 1 , y 2 ) be points in R 2 and T > 0 be a fixed real number. The conducting body D is a two-dimensional bounded domain in R 2 with smooth bounding surface Γ = ∂D, for example, C 2 -smooth is sufficient. The closure of the body D is D = D ∪ Γ. Composed with time we have the following cylinders D T = D × (0, T ] and Γ T = Γ × (0, T ], respectively. The closures of D T and Γ T are given by DT = D × [0, T ] and ΓT = Γ × [0, T ], respectively.

We are interested in constructing the solution u to the heat equation in the domain D T , supplied with initial and Dirichlet boundary conditions, that is u solves ∂u(x, t) ∂t -∆u(x, t) = 0, (x, t) ∈ D T ,

u(x, t) = h(x, t), (x, t) ∈ Γ T ,

u(x, 0) = u 0 (x), x ∈ D, (2)
where u 0 (x) and h(x, t) are sufficiently smooth functions. We point out that, in principle, the MFS that we propose and investigate can be applied to other boundary conditions, such as Neumann and mixed boundary conditions. To guarantee the existence and uniqueness of a solution to (1)-(3) we impose the following compatibility conditions:

u 0 (x) = h(x, 0) and ∂h ∂t (x, 0) = ∆u 0 (x), x ∈ Γ. ( 4
)
With these conditions the following uniqueness theorem holds, see, for example, [START_REF] Friedman | Partial Differential Equations of Parabolic Type[END_REF]. Theorem 2.1 Let u 0 (x) ∈ C 2 ( D) and h(x, t) ∈ C 1 ( ΓT ) satisfy the compatibility conditions (4). Then there exists a unique solution u ∈ C 2,1 ( DT ), to the equations (1)-( 3), which depends continuously on the data. Theorem 2.1 tells us, in particular, that the problem given by the equations (1)-( 3) is well-posed.

Denseness properties of linear combinations of fundamental solutions

The fundamental solution of (1) in two-dimensions is given by where H is the Heaviside function which is introduced in order to emphasize that the fundamental solution is zero for t ≤ τ . We shall investigate some properties of linear combinations of such functions for various source points y.

F (x, t; y, τ ) = H(t -τ ) 4π(t -τ ) e -|x-y| 2 4(t-τ ) , (5)
We begin by constructing a set of source points placed outside the region D. Let D E (E for enclosed) be an open domain, containing D, with bounding surface Γ E , where the distance between the points on the surfaces Γ and Γ E is greater than zero. Let {y j , τ m } j,m=1,2,... be a denumerable, everywhere dense set of points in Γ E × [-T, T ], (τ m = 0) and set

v (j) m (x, t) = F (x, t; y j , τ m ). ( 6
)
Figure 1 shows how the source points may be placed around a domain D, either using a symmetric shape or shapes which take the general shape of Γ obtained by dilatation. We now construct the following infinite series

u ∞ (x, t) = ∞ j=1 ∞ m=1 c (j) m v (j) m (x, t), (7)
where c

(j) m are set equal to zero except for a finite number of values. Note that, due to the Heaviside function in [START_REF] Chantasiriwan | The method of fundamental solutions for free surface Stefan problems[END_REF], we have u ∞ (x, t) = 0 for t ≤ τ = min m,j:|c (j) m | =0 τ m . Also note that, since F solves the heat equation, u ∞ also satisfies the heat equation in D T .

Denseness on the lateral surface

We prove the following denseness result on the lateral surface Γ × (-T, T ):

Theorem 3.1 The set of functions {v (j) m (x, t)} ∞ j,m=1
restricted on Γ×(-T, T ) form a linearly independent and dense set in L 2 (Γ × (-T, T )).

Proof A similar version of the proof of this theorem was given in one-dimension in [START_REF] Johansson | A method of fundamental solutions for transient heat conduction[END_REF] and in three-dimensions in [START_REF] Kupradze | A method for the approximate solution of limiting problems in mathematical physics[END_REF], and we follow those ideas here in the twodimensional case.

Linear independence: Assume that we do not have linear independence, then there exist positive integers N, m 0 , j 0 ∈ {1, . . . , N }, and a coefficient c

c (j) m v (j) m (x, t) = 0, (x, t) ∈ Γ × (-T, T ). ( 8
)
Define the function

U (x, t) = N j=1 N m=1 c (j) m v (j) m (x, t), (x, t) ∈ D × (-T, T ). ( 9
)
Then U satisfies the following equations:

∂U ∂t -∆U = 0, in D × (-T, T ), (10)
U (x, t) = 0, (x, t) ∈ Γ × (-T, T ), (11)
U (x, -T ) = 0. ( 12
)
We have obtained the above equations by observing that the fundamental solution satisfies [START_REF] Golberg | The method of fundamental solutions for potential, Helmholtz and diffusion problems[END_REF]; the Heaviside function makes the fundamental solution equal to zero in equation [START_REF] Gorzelanczyk | Some remarks concerning the shape of the source contour with application of the method of fundamental solutions to elastic torsion of prismatic rods[END_REF], and [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF] gives us [START_REF] Golberg | The method of fundamental solutions for time-dependent problems[END_REF]. By the uniqueness Theorem 2.1, the only solution to the equations ( 10)-( 12) is U (x, t) ≡ 0 for (x, t) ∈ D × (-T, T ). Because U is analytic in D E × (-T, T ), we also have U (x, t) = 0 for (x, t) ∈ D E × (-T, T ), see [START_REF] Saut | Unique continuation for some evolution equations[END_REF].

We now let the point (x, t) approach the point (y j0 , τ m0 ) ∈ Γ E × (-T, T ) such that the ratio |xy j0 | 2 4(t -τ m0 ) [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the L-curve[END_REF] remains bounded. Then the summand c 8) may be made as large as we wish, while the other terms in the series (8) remain bounded; this gives us a contradiction and thus, we have linear independence for the set of functions {v T,T )). Assume on the contrary that it is not a dense set. Then there exists an element f (x, t) in L 2 (Γ × (-T, T )), which we can assume is continuous, such that

(j0) m0 v (j0) m0 (x, t) in (
(j) m (x, t)} ∞ j,m=1 in L 2 (Γ × (-T, T )). Denseness: We next prove that the sequence {v (j) m (x, t)} ∞ j,m=1 is a dense set in L 2 (Γ × (-
T -T Γ v (j) m (x, t)f (x, t) dx dt = 0, j, m = 1, 2, . . . ( 14
)
To show that {v (j) m (x, t)} is dense we have to show that f (x, t) ≡ 0 in [START_REF] Hon | A discrepancy principle for the source points location in using the MFS for solving the BHCP[END_REF]. From definition [START_REF] Chen | Some comments on the ill-conditioning of the method of fundamental solutions[END_REF], equation ( 14) can be rewritten as where the Heaviside function in [START_REF] Chantasiriwan | The method of fundamental solutions for free surface Stefan problems[END_REF] has been used to reduce the range of integration with respect to the time variable. We introduce an equivalent form of the classical single-layer heat potential given by

T τm Γ F (x, t; y j , τ m )f (x, t) dx dt = 0, j, m = 1, 2, . . . , (15)
V (y, τ ) = T τ Γ F (x, t; y, τ )f (x, t) dx dt, y ∈ Γ. ( 16
)
It is well-known that V (y, τ ) is a smooth solution to the heat equation in the exterior of D × (-T, T ) and it cannot vanish on any surface in this exterior region without being identically zero. Thus, by the continuity of F and (15), we find that V (y, τ ) = 0 for (y, τ ) ∈ Γ E × (-T, T ), which is in the exterior of D × (-T, T ); we then conclude that V = 0 in the exterior of D × (-T, T ). Moreover, since V is continuous across Γ × [-T, T ] we also have V (y, τ ) = 0 on Γ × [-T, T ]. This implies that V = 0 also in D × (-T, T ) since V satisfies the heat equation in D × (-T, T ). Finally, using the jump relations for the normal derivative of

V on Γ × [-T, T ], see [9, p. 133], we get 1 2 f (x, t) ± ∂V (x, t) ∂ν = 0, (x, t) ∈ Γ × (-T, T ), ( 17
)
where ν represents the unit normal on the surface Γ × (-T, T ). Thus, f ≡ 0 and therefore, {v

(j) m (x, t)} ∞ j,m=1 is a dense set in L 2 (Γ × (-T, T )).

Denseness on the base surface

We now show that we also have denseness on the "base" surface D × {0}, where the initial condition is imposed in ( 1)-( 3).

Theorem 3.2

The set of functions {v

(j) m (x, 0)} ∞ j,m=1
, where v

m (x, t) is given by (6) with τ m < 0, form a linearly independent and dense set in L 2 (D).

Proof The method of proof is similar to that used in [START_REF] Johansson | A method of fundamental solutions for transient heat conduction[END_REF] in one-dimension, and we give it here, for completeness, in higher dimensions.

Linear independence: Assume that we do not have linear independence, then there exist positive integers N, m 0 , j 0 ∈ {1, . . . , N }, and a coefficient c

(j0) m0 = 0 such that N j=1 N m=1 c (j) m v (j) m (x, 0) = 0, x ∈ D. (18)
We shall use the corollary of Theorem 3 in [START_REF] Johnson | A Priori Estimates and Unique Continuation Theorems for Second Order Parabolic Equations[END_REF], which guarantees that if u(x, t) is a smooth solution of the heat equation in R 2 that is bounded by Be β|x| 2 and if u(x, 0) = ǫ then, |u(x, t)| ≤ ǫ, for 0 ≤ t ≤ T. From ( 9) and ( 18) we have that

U (x, 0) = 0, x ∈ D. ( 19
)
We also have that U satisfies [START_REF] Golberg | The method of fundamental solutions for potential, Helmholtz and diffusion problems[END_REF] and [START_REF] Gorzelanczyk | Some remarks concerning the shape of the source contour with application of the method of fundamental solutions to elastic torsion of prismatic rods[END_REF]. Now, since U (x 1 , x 2 , 0), where x = (x 1 , x 2 ), is a real analytic function in each of the variables x 1 and x 2 , we find that U (x, 0) = 0 for every x ∈ R 2 , see [20, p. 14]. Moreover, since each τ m < 0, U is continuous on R 2 × [0, T ] and is at least twice continuously differentiable in R 2 × [0, T ]. Furthermore, U also satisfies the heat equation ( 1), and the following inequality clearly holds

|U (x, t)| ≤ Be β|x| 2 , t ∈ [0, T ], (20)
for some positive constants B and β. Thus, from Theorem 3 of [START_REF] Johnson | A Priori Estimates and Unique Continuation Theorems for Second Order Parabolic Equations[END_REF], and its corollary, we conclude that U (x, t) = 0 in R 2 × [0, T ]. In particular, from [START_REF] Saut | Unique continuation for some evolution equations[END_REF], we may extend

U such that U (x, t) = 0 also in D E × [-T, T ].
We now let the point (x, t) approach the point (y j0 , τ m0 ) ∈ Γ E × [-T, 0] such that the ratio (13) remains bounded. Now, the summand c

(j0) m0 v (j0)
m0 (x, t) may be made as large as we wish, while the other terms in the series [START_REF] Karageorghis | Matrix decomposition MFS algorithms[END_REF] remain bounded; this gives a contradiction and we have linear independence for the set of functions {v

(j) m (x, 0)} in L 2 (D).
Denseness: We shall show that the set of functions {v (j) m (x, 0)}, where τ m < 0, is a dense set in L 2 (D). Assume that this is not a dense set, then there exists a

function f ∈ C 2 (D) such that D v (j) m (x, 0)f (x) dx = 0, j, m = 1, 2, . . . (21)
We let w be a weak solution of the heat equation [START_REF] Balakrishnan | The method of fundamental solutions for linear diffusionreaction equations[END_REF], see [START_REF] Evans | Partial Differential Equations[END_REF], with initial condition w(x, 0) = f (x) and boundary condition w(x, t) = 0 for (x, t) ∈ Γ T . We may transform (21) using Green's identities, see [START_REF] Lions | Non-homogeneous Boundary Value Problems and Applications[END_REF], into the following form

T 0 Γ v (j) m (x, t) ∂w(x, t) ∂ν dxdt = 0, j, m = 1, 2, . . .
where ν is the outward pointing unit normal to Γ. From Theorem 3.1 we know that {v (j) m (x, t)} restricted on Γ T is a dense set in L 2 (Γ T ), and we may conclude that the normal derivative of w is zero on Γ T . Therefore, both w and ∂w ∂ν are zero on Γ T . From [START_REF] Saut | Unique continuation for some evolution equations[END_REF], we conclude that w(x, t) = 0 for (x, t) ∈ DT ; hence f ≡ 0, and {v (j) m (x, 0)}, where τ m < 0, is a dense set in L 2 (D).

4.The MFS for the heat equation in two-dimensions

The denseness results proved in the previous section, Theorems 3.1 and 3.2, which involved linear combinations of the fundamental solution (5) of the heat equation ( 1), enable us to describe a method for approximating the solution to the problem (1)-( 3). We note that the MFS we propose may be applied to domains of general shape and size and the source points may also be placed arbitrarily, for example, placed symmetrically on a circular pseudo-boundary, or to match similarly the general shape of the domain D, see Figure 2 and [START_REF] Gorzelanczyk | Some remarks concerning the shape of the source contour with application of the method of fundamental solutions to elastic torsion of prismatic rods[END_REF]. The only restriction on the placement of the source points is that they are located on the boundary Γ E outside the domain D, and placed relatively close to D such that u has no singularity in D E × [0, T ]. Also, it might be more practical to take the sources on the interval (-ǫ, T ), where ǫ > 0 is small, instead of the full interval (-T, T ). However, some preliminary numerical investigations in [START_REF] Hon | A discrepancy principle for the source points location in using the MFS for solving the BHCP[END_REF] for the backward heat conduction problem showed that ǫ cannot be chosen too small if no loss in accuracy and stability is to be secured. We search for an approximation to the solution of equations ( 1)-( 3) in the following form:

u M,N (x, t) = 2M m=1 N j=1 c (j) m F (x, t; y j , τ m ), (x, t) ∈ DT . ( 22
)
For simplicity, let us describe the MFS in the case of circular domains. We shall also consider rectangular domains in the next section, see Example 3 and4.

We consider a two-dimensional circular domain D, with boundary Γ and radius r 0 > 0, centred at the origin, and let us place the source points (y j ) j=1,N on a circle r 0 + h, h > 0, also centred at the origin. The parameter h > 0 will be chosen such that the error at the lateral and base surfaces is minimized (viz maximum principle for the heat equation).

Take the time points (τ m ) m=1,2M (each in the interval (-T, T )) as given by

τ m = 2(m -M ) -1 2M T, m = 1, . . . , 2M,
and using polar coordinates, place the source points in space at

y j = (r 0 + h, θ j ) = r 0 + h, 2πj N , j = 1, . . . , N.
In polar coordinates equation ( 22) is now represented by

u M,N (r, θ, t) = 2M m=1 N j=1 c (j) m F (r, θ, t; r 0 + h, θ j , τ m ). (23)
We have located N × 2M source points in total outside the domain D and in time; we place the same number of collocation points in total on ΓT ∪ (D × {0}), the boundary in time and the domain at time t = 0. Of course, the location of source and collocation points may be chosen arbitrarily, here we choose points for ease of calculation. Let We have located N × (M + 1) collocation points on the boundary, the remaining N × (M -1) points will be located on D when t = 0. We consider M -1 circles of radius

t i = i M T, i = 0, . . . , M,
r l = l M 1 2 r 0 , l = 1, . . . , M -1,
where the square root has been introduced to spread the points out within the domain, and not to cluster them at the centre. We place N equally spaced points on each circle such that

(r l , θ k ) = r l , 2πk N , k = 1, . . . , N,
see Figure 3 for a detailed graphical representation of the position of the various source and collocation points given above. We now impose the boundary and initial conditions (2) and (3) so that we can determine the unknown coefficients c (j) m in [START_REF] Ramachandran | Method of fundamental solutions: singular value decomposition analysis[END_REF]. In polar coordinates we obtain the equations

u M,N (r 0 , θ k , t i ) = h(r 0 , θ k , t i ), (24)
u M,N (r l , θ k , 0) = u 0 (r l , θ k , 0), (25)
where k = 1, . . . , N, i = 0, . . . , M and l = 1, . . . , M -1.

The system of equations ( 24) and ( 25) contains N ×(M -1)+N ×(M +1) = 2M N equations and 2M N unknowns, therefore, we may obtain a unique solution. We can represent this system of equations as

Ac = g, ( 26
)
where c is the vector of unknowns c

(j)
m , g is the vector representing the values of the functions u 0 and h at the respective collocation points, and A is the matrix corresponding to the value of the fundamental solution at the points outlined above. For certain boundary collocation and source points it might be possible to use the properties of circulant matrices to develop a matrix decomposition algorithm [START_REF] Karageorghis | Matrix decomposition MFS algorithms[END_REF] for the solution of the system of equations [START_REF] Young | Time-dependent fundamental solutions for homogeneous diffusion problems[END_REF] and thus, substantially reduce the computational cost. This could be considered in a future investigation.

To solve this system it is possible to solve it directly by Gaussian elimination. However, it is well-known that employing the MFS can yield matrices with large condition numbers [START_REF] Chen | Some comments on the ill-conditioning of the method of fundamental solutions[END_REF][START_REF] Ramachandran | Method of fundamental solutions: singular value decomposition analysis[END_REF], as h increases. Then, in such a situation instead of [START_REF] Young | Time-dependent fundamental solutions for homogeneous diffusion problems[END_REF] it may be necessary to consider the Tikhonov regularization

(A tr A + λI)c = A tr g, (27)
where the superscript tr denotes the transpose of a matrix and I is the identity matrix. The linear system of equations ( 27) is solved using a Gaussian elimination method (employed backslash "\" command in MATLAB). In (27), λ > 0 is a small regularization parameter (usually in the interval [10 -1 , 10 -16 ]) can be chosen by trial and error, namely start with a large value of λ, say λ = 10 -1 , and then decrease it gradually as λ = 10 -2 , 10 -3 , . . . until an oscillatory solution starts to develop. To choose λ one can also use the L-curve criterion of [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the L-curve[END_REF]. In a future study, it would be interesting to investigate the dependence of the solution on λ and h.

Numerical results

In [START_REF] Johansson | A method of fundamental solutions for transient heat conduction[END_REF] it was shown that the direct MFS approximation applied to the onedimensional heat equation with source points located outside the domain and in time is accurate. Below, we present numerical results for approximations in twodimensional domains, such as circular and square domains. In order to assess the accuracy of the numerical MFS solutions we compare them with the available exact solutions for various benchmark test examples. Numerical results are presented for N = 20 and M = 30 points, which were found sufficiently large to ensure that any further increase in these numbers did not significantly improve the accuracy of the numerical solution without affecting its stability.

Example 1

Let D = {x : |x| 2 < 1}, D T = {(x, t) : |x| 2 < 1, t ∈ (0, 1]}, and Γ T = {(x, t) : |x| 2 = 1, t ∈ (0, 1]}. We solve the following problem, using the direct MFS laid out in the previous section,

∂u(x, t) ∂t -∆u(x, t) = 0, (x, t) ∈ D T , (28)
u(x, t) = 4t + 1, (x, t) ∈ Γ, ( 29
) u(x, 0) = |x| 2 , x ∈ D. (30)
The exact solution of problem (28)-(30) is u(x, t) = 4t + |x| 2 . The source points are placed on a circle with radius 1 + h. The value of h > 0 will be chosen appropriately. However, the accuracy of the approximation appears to decrease when h < 0.25 or h > 4. In Figure 4 the exact solution and the MFS approximations are plotted in one-dimension, x = (x 1 , 0), for times t ∈ {0.2, 0.8} with λ = 10 -8 in the Tikhonov regularization. As can be seen in Figure 4 the approximation obtained is accurate, and was stable even when considering other small positive values of λ, however, for stable results we could not take λ = 0 due to the round-off precision errors.

-1 -0.5 0 0.5 1 0.5 Finally, we consider a three-dimensional plot of the exact solution u(x 1 , x 2 , 0.8) in Figure 6(a), and the MFS approximation u M,N in Figure 6(b) obtained with h = 1. Figure 6(c) shows the graph of the absolute error, and we note that the approximation is very accurate with a maximum absolute error of O(10 -5 ).

Example 2

In this example we choose the same D, D T and Γ as in Example 1, but instead of u(x, t) = 4t + |x| 2 , we consider the exact solution of the equations ( 1)-(3) given by u(x, t) = e x1+x2 cos(

x 1 + x 2 + 4t), ( 31
)
where the boundary and initial equations ( 2) and ( 3) have been obtained from (31). Note that this function is not constant on circles centred at the origin, and not symmetric, thus being different in character compared with the solution in Example 1.

Again, the source points will be placed on a circle with radius 1 + h, with final time point T = 1. As an alternative to the strategy for choosing λ described at the end of Section 4.1, we employ the L-curve criterion [START_REF] Hansen | Analysis of discrete ill-posed problems by means of the L-curve[END_REF]. In Figure 7 we present a plot of the L-curve for Example 2, where the residual is plotted against the 2-norm of the solution c. We choose the regularization parameter λ, which corresponds to the corner of the "L" in Figure 7; namely, in this example we take λ = 10 -8 . In Figure 8 the exact solution and the direct MFS approximations for λ = 10 -8 are plotted in one-dimension, x = (x 1 , 0), for times t = {0.2, 0.8}, with h = 1. It is clear from Figure 8 that the approximation is very accurate. Also, varying λ between [10 -12 , 10 -4 ] only slightly changes the accuracy of the approximation. Finally, we consider a three-dimensional plot of the exact solution u(x 1 , x 2 , 0.8) The plots obtained in Examples 1 and 2 show that the MFS approximation is accurate in circular domains, with errors usually in the interval [10 -4 , 10 -2 ], for a wide range of parameters h and λ.

1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -3 -2 -1 0 1 x 1 x 2 u(x 1 , x 2 , 0.2) (a) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -3 -2 -1 0 1 x 1 x 2 u(x 1 , x 2 , 0.2) (b) -1 0 1 -1 0 1 0 0.5 1 1.5 x 10 -4 x 1 x 2 u error (x 1 , x 2 , 0 .2) (c)

Example 3

In the next two examples we consider square domains with edge length L. In Examples 3 and 4 source points will be placed on both squares, see Figure 10, as well as circles; we vary the shapes where we place the source points to highlight that the placement of the source points do not need to follow the shape of Γ. The numerical implementation is the same as in the previous examples, including the placement of the sources and boundary collocation points and their numbers N = 20, M = 30. The following problem was considered in [START_REF] Golberg | The method of fundamental solutions for time-dependent problems[END_REF]. We take D = (-0.2, 0.2) × (-0.2, 0.2) to be a square of edge length L = 0.4, take T = 0.9, and solve

x 1 x 2 t T -T Γ E Γ D L L + h
∂u(x, t) ∂t -α∆u(x, t) = 0, (x, t) ∈ D T = D × (0, 0.9], (32)
u(x, t) = 0, (x, t) ∈ Γ T = ∂D × (0, 0.9], (33) Here, we have thermal diffusivity α = 1000/5.8, and the fundamental solution of equation ( 32) is now given by

u(x, 0) = 1, x ∈ D. (34)
F (x, t; y, τ ) = H(t -τ ) 4πα(t -τ ) e -|x-y| 2 4α(t-τ ) . (35)
We note that in this example the compatibility conditions (4) are violated. The exact solution to the problem (32)-( 34) is given by, see [START_REF] Carslaw | Conduction of Heat in Solids[END_REF],

u(x 1 , x 2 , t) = 16 π 2 ∞ n=0 (-1) n 2n + 1 e -α(2n+1) 2 tπ 2 /(4(L/2) 2 ) cos (2n + 1)πx 1 2(L/2) × ∞ m=0 (-1) m 2m + 1 e -α(2m+1) 2 tπ 2 /(4(L/2) 2 ) cos (2m + 1)πx 2 2(L/2) . (36)
However, when we plot the exact (using 100 terms in the series expansion (36)) and the approximate solutions at the point (x 1 , x 2 ) = (0, 0) for time t ∈ [0, 0.9], with source points placed on a square with edge length 0.4 + h, h = 5, λ = 10 -8 , we get a large discrepancy, see Figure 11. Changing the parameter λ in the Tikhonov regularization (27) did not seem to improve the approximation, however, we observe that the exact solution (36) decays very rapidly due to the exponential terms. This means that we should consider a much smaller time interval. Figure 12 shows the exact and the approximate solutions, with final time point T = 0.0006, plotted over t ∈ [0, 0.0004], where λ = 10 -6 and source points have now been placed on a circle of radius h = 0.84.

These figures show that the choice of the final time T , in particular when considering fast decaying functions, is also important when implementing the MFS. Time-marching methods, [START_REF] Valtchev | A time-marching MFS scheme for heat conduction problems[END_REF], could then perhaps be used to extend our approximation to larger time intervals.

Example 4

We consider D = (0, 1) × (0, 1) and D T = D × (0, 3], and we wish to solve the following problem:

∂u(x, t) ∂t -∆u(x, t) = 0, (x, t) ∈ D T , (37)
u(x 1 , 0, t) = u(x 1 , 1, t) = √ 2e -π 2 t/4 cos πx 1 2 - π 4 + 1 √ 2 , x 1 ∈ (0, 1), t ∈ (0, 3], (38)
u(0, x 2 , t) = √ 2e -π 2 t/4 cos πx 2 2 - π 4 + 1 √ 2 , x 2 ∈ (0, 1), t ∈ (0, 3], ( 39
)
∂u ∂x 1 (0, x 2 , t) = π 2 e -π 2 t/4 , x 2 ∈ (0, 1), t ∈ (0, 3], (40) u
(x 1 , x 2 , 0) = √ 2 cos πx 1 2 - π 4 + cos πx 2 2 - π 4 , (x 1 , x 2 ) ∈ D. (41)
This is an inverse problem with missing boundary data at x 1 = 1, which we wish to determine using the Cauchy boundary data over-specification at x 1 = 0. The exact solution of problem (37)-(41) is

u(x 1 , x 2 , t) = √ 2e -π 2 t/4 cos πx 1 2 - π 4 + cos πx 2 2 - π 4 , (x, t) ∈ D T . (42)
In this example, we shall show results when the source points are placed on squares, as well as when we place them on circles, to show that placement of the sources does not need to follow the shape of the solution domain. When the source points are placed on a square they will be located at (-h/2, 1 + h/2) × (-h/2, 1 + h/2), whilst the source points placed on a circle will have radius h with centre at (0.5, 0.5), see Figure 13 Figure 15 shows the exact solution u(x 1 , x 2 , 0.5) and the MFS approximation u M,N . Note that for the approximation u M,N in Figure 15(c) we have instead placed the source points on a circle and there is still good agreement with the exact solution. In Figure 16 we present plots of the absolute error at time t = 2.5 for two different values of h. From this figure it can be seen that the error increases when source points have been placed too close to the boundary. Finally, in Figure 17, random noise simulating measurement errors, have been added to the Dirichlet boundary data (39) as follows:

u δ (0, x 2 , t) = u(0, x 2 , t) + N (0, σ 2 ),
where N (0, σ 2 ) represents the normal distribution with mean zero and standard |u(0, x 2 , t)|, and δ is the relative noise level. A set of ten noisy random data functions {u δ k (0, x 2 , t)} k=1,10 was generated, and the source points in the MFS have been placed on a circle of radius h, centred at (1/2, 1/2). Figure 17(a) presents a plot of the exact solution u(1, x 2 , 0.5), and the best ( * ) and the least accurate (•) MFS approximations from these ten data sets, obtained with δ = 3% noise, h = 3, T = 3 and λ = 10 -4 . In Figure 17(b) we present a three-dimensional plot of the absolute error when δ = 3%, h = 3, λ = 10 -4 and t = 0.5. Adding more noise such as δ = 5% did not significantly change the stability of the numerical results provided that regularization is applied appropriately. As expected, the accuracy decreases when the noise level increases and the regularizing parameter usually has to take a larger value. Noise can also be added to the other data functions, such as the Neumann data (40), and the same stable and accurate numerical results are expected. Thus, for this inverse problem, the regularized MFS is a stable approximation with respect to noisy data.

Conclusions

We have investigated the application of the MFS for linear heat conduction in twodimensional conducting bodies with Dirichlet boundary conditions. The solution is sought in the form of a linear combination of fundamental solutions of the heat equation with the source points placed outside the body and in time as well. Theoretical properties, such as denseness of the approximation on the boundary of the body, have been obtained in the space of square integrable functions. Moreover, numerical results are presented for both circular and rectangular configurations showing that accurate approximations can be obtained when the source points are placed on circles or rectangles at an appropriate distance from the solution domain. The numerical results show that accurate approximations can be obtained at small computational cost, and that the method is not too sensitive with respect to the placement of the sources. Numerical results have also been included for an inverse problem, where overspecified data was given on a part of the boundary of a body and the solution was to be reconstructed on the remaining part. Good agreement with the exact available analytical solution was also obtained in this case, showing that boundary conditions other than Dirichlet can be handled. Extensions to the three-dimensional time-dependent heat equation are straightforward by accordingly modifying the fundamental solution [START_REF] Chantasiriwan | The method of fundamental solutions for free surface Stefan problems[END_REF].

Figure 1 .

 1 Figure 1. MFS for two-dimensional heat conduction, source points located outside of the spatial domain D and in time [-T, T ].

  Tomas Johansson et al.

Figure 2 .

 2 Figure 2. Arbitrary domain with varying source point locations restricted to any contour Γ E embracing the given solution domain D.4.1 A direct MFS for the two-dimensional heat equation

Figure 3 .

 3 Figure 3. Location of source and collocation points when D is a circular domain.

Figure 4 .

 4 Figure 4. The exact solution (-) and the approximate values ( * ) for u(x 1 , 0, 0.2) and u(x 1 , 0, 0.8), for Example 1.

Figure 5 Figure 5 .

 55 Figure5contains plots of the exact solution and the direct MFS approximations for h ∈ {0.5, 4} and λ = 10 -8 . From this figure it can be seen that the numerical results obtained with h = 0.5 are slightly more accurate than those obtained with h = 4.

Figure 6 .

 6 Figure 6. Three-dimensional plot of: (a) The exact solution u(x 1 , x 2 , 0.8), (b) the approximate solution u M,N obtained with λ = 10 -8 and h = 1, and (c) the absolute error, for Example 1.

Ac -g 2 c 2 λ = 10 - 1 λ

 22101 = 10 -5 λ = 10 -8 λ = 10 -10 λ = 10 -12

Figure 7 .

 7 Figure 7. Plot of the L-curve, for Example 2.

1 u

 1 (x 1 , 0, 0.2) u(x 1 , 0, 0.8)

Figure 8 .

 8 Figure 8. The exact (-) and the approximate values ( * ) for u(x 1 , 0, 0.2) and u(x 1 , 0, 0.8), for Example 2.

  a), the MFS approximation u M,N in Figure9(b), and the absolute error in Figure9(c), obtained with h = 1.

-

Figure 9 .

 9 Figure 9. Three-dimensional plot of: (a) The exact solution u(x 1 , x 2 , 0.2), (b) the approximate solution u M,N obtained with λ = 10 -8 and h = 1, and (c) the absolute error, for Example 2.

Figure 10 .

 10 Figure 10. MFS in two-dimensions for square domains.

Figure 11 .

 11 Figure 11. The exact solution u(0, 0, t) and the approximation u M,N , as functions of time t ∈ [0, 0.9], for Example 3.

Figure 12 .

 12 Figure 12. The exact solution u(0, 0, t) and the approximation u M,N with T = 0.0006, for Example 3.

  .

Figure 13 .

 13 Figure 13. Examples of source point location for the unit square solution domain of Example 4.

Figures 14 (Figure 14 . 1 ( 1 ,

 141411 Figures14(a) and14(b) show the exact solution u(1, x 2 , 1.5) and its normal derivative ∂u ∂x1 (1, x 2 , 1.5), respectively, in comparison with the approximate solutions, obtained with h = 3 and λ = 10 -8 .

Figure 15 .Figure 16 .

 1516 Figure 15. (a) The exact solution u(x 1 , x 2 , 0.5), (b) the MFS approximation u M,N using source points placed on a square, h = 3, and (c) the MFS approximation u M,N using source points placed on a circle with radius h = 3 and centre (0.5, 0.5), for Example 4.

Figure 17 .

 17 Figure 17. (a) The exact solution u(1, x 2 , 0.5) (-) and the best ( * ) and the least accurate (•) MFS approximations from ten different sets of noisy data with noise level δ = 3% and (b) the absolute error |u(x 1 , x 2 , 0.5)u M,N (x 1 , x 2 , 0.5)| when δ = 3%, for one of the noisy data sets, for Example 4.

Acknowledgements

The authors would like to thank the referees for their suggestions. The third author acknowledges financial support from EPSRC in the form of a Doctoral Training Grant (DTG).


Sandra Engel
email: sandra.engel@uos.de

Anna Alemany

Nuria Forns

Philipp Maass
email: philipp.maass@uos.de

Felix Ritort
email: fritort@gmail.com

Folding and unfolding of a triple-branch DNA molecule with four conformational states

Keywords: biophysics, statistical physics, biomolecules

Single-molecule experiments provide new insights into biological processes hitherto not accessible by measurements performed on bulk systems. We report on a study of the kinetics of a triple-branch DNA molecule with four conformational states by pulling experiments with optical tweezers and theoretical modelling. Three distinct force rips associated with different transitions between the conformational states are observed in the folding and unfolding trajectories. By applying transition rate theory to a free energy model of the molecule, probability distributions for the first rupture forces of the different transitions are calculated. Good agreement of the theoretical predictions with the experimental findings is achieved. Furthermore, due to our specific design of the molecule, we found a useful method to identify permanently frayed molecules by estimating the number of opened basepairs from the measured force jump values.

Introduction

In recent years, single-molecule experiments became of great importance in biophysical research since progress in nano-and microscale manufacturing technologies facilitated the design of scientific instruments with sufficient sensitivity and precision to enable the controlled manipulation of individual molecules (for reviews, see, for example, [1,2]). In contrast to the traditionally used bulk assays, where individual biomolecular dynamics can get masked, single-molecule experiments provide new insights into the thermodynamics and kinetics of biophysical and biochemical processes hitherto not accessible. They complement standard spectroscopy and microscopy methods used in molecular biology and biochemistry and hence have to be regarded as an important source of additional information helping in the interpretation of biomolecular processes. Furthermore, single-molecule experiments permit the measurement of small energies and the detection of large fluctuations.

The manipulation of single molecules offers a powerful new tool in molecular and cellular biophysics allowing for the exploration of processes occuring inside the cell at an unprecedented level. To instance just a few of the recently investigated biochemical processes: the transport of matter through pores or channels [3][4][5], interactions between DNA and proteins [6] or DNA and RNA [7], the motion of single-molecular motors [8][9][10], DNA transcription and replication [11,12], virus infection [13,14], DNA condensation [15] and ATP generation [16]. In addition, the structure of biological networks [17] and the viscoelastic and rheological properties of the DNA [18][19][20][21] have been studied. An important class of single-molecule experiments are performed with optical tweezers. By means of an optical trap generated by a focused laser beam, this useful technique renders it possible to exert forces on micron sized objects, achieving sub-piconewton and sub-nanometer resolution in force and extension, respectively. Accordingly, one can study force-induced folding-unfolding dynamics and in this way get insight into corresponding processes in the cell and typical bond forces. Of particular interest is the unfolding of DNA molecules, where the hydrogen bonds between the complementary base pairs (bps) are disrupted. This so-called unzipping is connected to the DNA replication mechanism.

An interesting field of biophysical studies is the investigation of junctions in molecules since they present manifold ways to interact with other substances, for instance cations. Three-way junctions are especially interesting because metal ions such as magnesium can bind to them and alter the tertiary structure. Here a first step of such a study is presented where we investigate a molecule with a three-way junction alone, without cation binding.

Many of the single-molecule experiments so far focused on molecules with a relatively simple free energy lanscape (FEL) exhibiting just two states, a folded and an unfolded one, or including an additional misfolding state, leading to different kinetic pathways. In this work we will consider a richer situation, where metastable states as intermediates occur during the folding-unfolding route. In this context, we will address the following key questions:

(i) Can a corresponding molecule with such intermediate states be designed on the basis of a suitable model for a FEL? (ii) Is it possible to observe the intermediate states by perfoming pulling experiments with optical tweezers? (iii) Can phenomenological Bell-Evans kinetic models be applied to describe the folding-unfolding processes including intermediate states? In particular, when validating the kinetic theory against the experimental results, how do the first rupture force distributions compare with the ones predicted by the theory?

A further important aspect that we looked at in some detail is the heterogenity of molecular folding-unfolding behaviour that we observed in the experiments. By measuring force-distance curves (FDCs) of several molecules we classified them into different reproducible patterns. This leads to a useful method to identify irreversible molecular fraying, a phenomenon which is often observed in single-molecule studies.

Description of the experiments

Based on Mfold folding predictions [START_REF] Santalucia | Proc. Natl. Acad. Sci[END_REF][START_REF] Zuker | [END_REF] and taking FEL considerations into account (see sec. 5), we designed and synthesised a DNA molecule which is composed of three parts and hence referred to as triple-branch molecule. It consists of a stem as introduced in [24] with 21 bps and two nearly identical hairpin branches which are formed of 16 bps and a loop with four bases, thus comprising a total number of 114 bases, see fig. 1. To avoid misfolding, the second hairpin branch differs at two positions from the first one.

The triple-branch molecule is inserted between two identical short double-stranded DNA (dsDNA) handles of 29 bps each [25], leading to a total number of 172 bases, corresponding to a total contour length in the unfolded state of about 100 nm. Each of these polymer spacers is chemically linked to a bead. One of the latter is retained with the help of a pipette via air suction, the other one is optically trapped in a laser focus [26]. On the 5 end of the DNA molecule biotin is attached to enable a connection with a streptavidincoated bead (SA bead) whose diameter is 1.8 µm. Biotin is a vitamin which establishes a strong linkage to the proteins avidin and streptavidin. The 3 end is modified with the antigen digoxigenin, able to interact with an antidigoxigenin-coated bead (AD bead). The latter has a diameter of 3.0 µm. Figure 2 shows the different components of the molecu- lar construct, that is to say the triple-branch DNA molecule, the handles and the beads, captured in optical trap and micropipette, respectively. Note that it is not a true-to-scale representation.

The pulling experiments are carried out with a miniaturised dual-beam laser optical tweezers apparatus [27] at room temperature ( 25 °C) and at salt concentration of 1 M NaCl aqueous buffer with neutral pH (7.5) stabilised by Tris HCl and 1 M EDTA. The dual-beam optical tweezers collect data at 4 kHz and can operate with a feedback rate of 1 kHz. Spatial resolution constitutes 0.5 nm with a maximal distance range of 10 µm. Forces up to 100 pN can be achieved, whereas the force resolution is 0.05 pN.

Optical tweezer pulling experiments permit the measurement of the force f as well as the total distance X between the centre of the optical trap and the tip of the micropipette, see fig. 2. In the experiment we vary the trap-pipette distance X(t) with a constant speed v = dX/dt in the range of 45 nm/s to 200 nm/s, which corresponds to a constant average loading rate r of 3.0 pN/s and 13.4 pN/s in between rip events, respectively. The experiment consists of loading cycles which in turn are divided into an unfolding part (during "pulling") and a folding part (during "pushing"). The loading cycles are repeated as long as the tether connection is unbroken. Otherwise a new connection has to be established, possibly a new molecule must be searched and linked to a new bead. In sec. 4 we will discuss different patterns found in the measured curves and present a detailed analysis of two representative molecules. We chose them among 7 molecules exhibiting the first and 5 molecules featuring the second pattern. For each molecule we recorded, on average, approximately 50 cycles. The pulling speeds, ranging from 45 to 200 nm/s, influence the experimental results only weekly due to a logarithmic dependence of the first rupture force with the speed. We found compatible data for sets of similar molecules. In the theoretical analysis of the data in sec. 5 we concentrate on the largest set of 82 loading cycles for the molecule at a speed of 200 nm/s. The study of the above mentioned second molecule comprises 55 cycles.

Analysis of unfolding and folding trajectories

Based on the design of the triple-branch molecule, we have to distinguish between four conformational states (see fig. 3):

(1) a completely folded molecule.

(2) a completely unfolded stem, whereas the hairpin branches are still folded.

(3) stem and either hairpin branch 1 or 2 are completely unfolded. (4) a completely unfolded molecule. Figure 4(b) displays a typical unfolding and folding trajectory in form of a FDC which in fact records the evolution of the force as a function of time and relative trap position X1 . With rising X the force f first increases almost linearly according to an elastic response of the DNA handles, which consist of dsDNA and are stable over the whole range of forces where the unfolding / folding of the triple-branch molecule takes place. The overstretching transition of the linkers, typically at 65 pN, lies much above the forces we explore. At a first rupture force f 1 a sudden decrease ("jump") ∆ f 1 occurs which is caused by the unfolding of the stem. This unfolding goes along with an abrupt change in the length of the molecule when single-stranded DNA (ssDNA) is released. As a consequence, the bead in the optical trap moves towards the centre of the trap, visible as the force drops. Following the jump ∆ f 1 , there is again a linear increase up to the next force rip at a first rupture force f 2 , where one of the hairpin branches unfolds, which in turn leads to the force jump ∆ f 2 . Eventually, the second hairpin branch unzips at a first rupture force f 3 with a jump ∆ f 3 . The linear regime following this last force rip corresponds to the stretching of the whole molecular construct including handles and the already unfolded triple-branch molecule.

Upon decreasing X from the completely unfolded state 4, the force first follows closely the corresponding unfolding part of the trajectory. However, a backward transition does not occur at ( f 3 -∆ f 3 ) 16.2 pN, but at a considerably lower value of (14.0 ± 0.8) pN, cp. fig. 4(a), hence manifesting a hysteresis effect. Moreover, an investigation of a larger number of folding trajectories reveals that the bases of the hairpin branches do not always pair conjointly in well-defined events during a short time interval. In contrast, the corresponding force rips during unfolding indicate a cooperative behaviour of the biomolecule, where the breakage of all hydrogen bonds stabilising the DNA structure happens almost simultaneously. The second transition takes place at (13.3 ± 0.8) pN instead of ( f 2 -∆ f 2 ) 15.5 pN. After both hairpin branches refolded, one can identify another sharp transition to the folded state 1 around (11.2 ± 0.8) pN, which is again considerably lower than ( f 1 -∆ f 1 ) 13.7 pN.

During unfolding, one can assume that the breakage of hydrogen bonds follows the sequence of the molecular construct. This makes it useful to introduce the number n of broken bonds as state variable and to calculate a FEL as function of this variable (see sec. 5). The refolding of the hairpin branches in the folding trajectories exhibit less sharp transitions, see fig. 4(b). During folding, in particular at the beginning in the unfolded state, a huge number of secondary structures can be found which implies that the kinetic pathways are less predefined and accordingly, the transitions get smeared out. With respect to a theoretical treatment, moreover, a description in terms of the simple state variable n becomes unlikely to be sufficient. In a refined analysis, many more configurations should have to be included as relevant states in a coarse-grained description [28]. Such refined analysis, however, goes beyond the scope of this work amd we therefore concentrate on the unfolding process in the following.

There are plenty of possible ways to analyse the FDCs in order to find out the first rupture forces and the force jumps. In our procedure we arranged the normalised data, i.e. the relative distance X and the force f , in windows of a certain size of data points. For all consecutive windows we then calculated the slope of the considered data points, the span, i.e. the maximum distance between the lowest and the highest force value, and the mean of the relative distance X j as a moving average. Transitions between the conformational states take place where the slope is minimal and the span maximal under the condition that an appropriate number of contiguous windows is connected. Having found the X j of the three force rips, slope and axis intercept are calculated by linear regression for each conformational state. One can now easily calculate the first rupture forces as the intersection points with the four fitted lines and extract the force jump values, as exemplified in fig. 5 for both molecules whose unfolding trajectories were depicted in figs. 4(a) and 4(b), respectively. Since the data acquisition rate is constant, at higher pulling speeds less data points are collected and therefore the values of the analysis are broader distributed. The three first rupture forces f i as well as the jumps ∆ f i are subject to stochastic fluctuations, as can be seen in fig. 4(a) where we show four unfolding trajectories belonging to different pulling cycles of the same molecule in the inset. An analysis revealed that the fluctuations of the force jumps ∆ f i are about ten times smaller than the fluctuations of the f i . Accordingly, in sec. 5, we will disregard the fluctuations in the ∆ f i and use only their averages ∆ f i that will be discussed in more detail in the following section. Applying the above mentioned procedure to analyse the experimental data of several molecules, we found two predominating patterns in the unfolding trajectories which are reflected in the distributions of the first rupture forces f i as follows. In the first pattern, see fig. 6(a), these distributions have a similar shape for all three force rips. The histograms indicate the existence of one maximum slightly below 17 pN. The mean values for the three first rupture forces are f 1 = (17.0 ± 0.9) pN, f 2 = (16.4 ± 0.6) pN and f 3 = (16.8 ± 0.7) pN. In contrast, in a second pattern we detected a strikingly lower value for the first rupture force of the first rip f 1 = (14.6 ± 0.8) pN, as depicted in fig. 6(b), whereas the other two rip forces lie basically in the same range of 16 to 18 pN. As in the former case, the second rip tends to have a slightly smaller first rupture force, f 2 = (16.5 ± 0.8) pN, than the third rip, f 3 = (17.2 ± 0.5) pN. The small f 1 observed in the second molecule suggests the occurrence of permanent molecular fraying. Obviously less force is needed to unfold the stem than typically, cp. fig. 4(a) with 4(b), since some bps at its basis are partly or completely melted. In other words, during folding, this molecule does not reach an entirely folded state but some bps of the stem next to the handles remain irreversibly and permanently open. This phenomenon has been observed previously in several pulling experiments [27,29]. A possible reason for irreversible fraying is the formation of reactive oxidative species due to the impact of the laser light of the optical trap leading to a degradation of the DNA bases [30].

The so generated singlet oxygens are known to oxidise certain nucleic acids, such as guanine and thymine, irreversibly. This could explain our observation that once a molecule shows fraying it does not change back again to normal behaviour. Due to the fact that we work with polystyrene microspheres which are more prone to photodamage than the DNA bases themselves, their wide ranging interaction with the bases might be reduced replacing polystyrene by silica beads which exhibit considerably minor irreversible oxidative damage. It would be very interesting to carry out such experiments. The insets in figs. 6(a) and 6(b) depict the corresponding force jump distributions of both molecules. While the first molecule possesses a large first force jump value of ∆ f 1 = (1.3 ± 0.2) pN and two smaller force jumps at the second and third rip of ∆ f 2 = (0.9 ± 0.2) pN and ∆ f 3 = (1.0 ± 0.2) pN, the frayed molecule features three force jumps at approximately the same value, i.e. ∆ f 1 = (0.93 ± 0.05) pN, ∆ f 2 = (0.95 ± 0.06) pN and ∆ f 3 = (0.97 ± 0.04) pN, respectively. This illustrates clearly the influence of irreversible fraying in the latter case since, roughly estimated, the same number of bps is expected to open in all three rips which should not be the case in an entirely folded molecule.

During a force rip, the relative distance X is constant and thus ∆X = 0. Therefore the change in the relative molecular extension,

∆x m = ∆X -∆ f /k eff = ∆x(n, f ), (1)
is only related to the force jump ∆ f and the combined stiffness of bead and handles k eff , given by 1

/k eff = 1/k b + 1/k h ,
where k b is the trap stiffness and k h the rigidity of the handles, respectively. With the help of a linear least squares fit, the effective stiffness k eff of the molecular construct is extracted from the average slope of the FDCs and amounts to (0.067 ± 0.003) pN/nm. Note that, as expected, k eff is smaller than k b ( 0.08 pN/nm). For a certain force value f and assuming an elastic model for the released ssDNA, the number n of opened bps is related univocally to the equilibrium end-to-end distance of the DNA molecule x(n, f ), whereas its change ∆x(n, f ), in turn, equals ∆x m . Using this relation, it is now possible to estimate the change in the number ∆n i of bps which are opened sequentially during each force rip. Considering only the configurations of the four conformational states, the equilibrium end-to-end distance x(n, f ) can be decomposed into two parts, cp. fig. 3. The first part, the elongation u l ( f ) of the mean end-to-end distance of the ssDNA along the force direction, accounts for the ideal elastic response of the ssDNA, where l is the contour length. The second part contains the contribution of the diameter of stem and hairpin branches, respectively. Acccordingly,

x(n, f ) = u l ( f ) +          0, n =21 + 16 + 16 bps d 0 , n = 0 and 21 + 16 bps d 0 , n =21 bps , (2)
where d 0 2 nm, in accordance with the diameter of the B-DNA helix. The exact value of the diameter contribution of both hairpin branches d 0 , when the stem is unfolded, depends on the orientation of the branches. We set d 0 = 2d 0 as a working value 1 .

Regarding the contour length l, which depends, amongst others (cp. sec. 5), on the number n of opened bps, and considering again solely the configurations of the four conformational states, one gets

l = 2nd +          0, n = 0
where the interphosphate distance d is taken to be 0.59 nm/base and n loop = 4 is the number of bases per loop. For every force rip f i the corresponding number of opened bps n i is calculated separately by considering the differences of l between the states. Different types of models can be used to calculate u l ( f ). Prominent examples borrowed from polymer physics are the freely jointed chain (FJC) and the worm-like chain (WLC) model. According to ref. [18], the FJC model includes an extra term, leading to the expression , respectively. In the WLC model [31], the force f (u l ), due to an elongation u l , is given by

u l ( f ) = l 1 + f Y coth b f k B T - k B T b f . ( 4
f (u l ) = k B T P 1 4 (1 -u l /l) 2 - 1 4 + u l l , (5)
and to obtain u l ( f ), this equation has to be inverted. Based on the WLC model, we tested the influence of the persistence length P in a typical range of 1.0 to 1.5 nm [20].

From the data shown in table 1 it is apparent that, depending on the model and parameters, the results for the estimated change in the number of opened bps vary in an acceptable range when the errors are taken into account. In addition, one can see that it is not evident which model and parameters should be considered as the best ones. Good results are found for the FJC model using recent values of [27] with P = 1.3 nm [32]. We chose to work with the FJC model with the parameters of [27] in sec. 5.

The second molecule indeed reveals a considerably smaller ∆n 1 , depending on the model around 13 or 14 bps, so that 7 or 8 bps are not closed after the folding process is completed. Performing single-molecule experiments without knowing the exact influence of permanently frayed bps can lead to misinterpreted results. Checking the appropriate parameters for the polymer models with the help of the change in the number of opened bps of the hairpin branches, one can estimate the number of irreversibly frayed bps at the basis of the stem1 .

Theory for the unfolding kinetics

The kinetics of the unfolding process can be described on a coarse-grained level based on a Gibbs free energy G(n, f ) as a function of the number n of sequentially opened bps for an applied force f . For small forces, including f = 0, the FEL is expected to have a shape as displayed in fig. 7(a). In general, G(n, 0) increases monotonously with n. However, local minima occur at the metastable states 2, 3 and 4 because there is an increase of entropy associated with the release of additional degrees of freedom when the stem-hairpin-junction and end-loops of the hairpin branches are opened. With rising force the FEL is expected to get tilted, so that the energies of the metastable states are lowered. With a knowledge of G(n, f ) we can apply standard transition rate theory and write for the transition rate from state i to i + 1

Γ i,i+1 ( f ) = γ 0 i γ i,i+1 ( f ) , (6)
where γ 0 i is an attempt rate and γ i,i+1 ( f ) is the Boltzmann factor corresponding to the activation barrier ∆G i,i+1 ( f ) that has to be surmounted

γ i,i+1 ( f ) = exp - ∆G i,i+1 ( f ) k B T . ( 7
)
When considering only sequential configurations in the evaluation of the FEL, different structures compatible with a given n can occur once the stem is completely unfolded. These refer to different possibilities of breaking the bps in the hairpin branches 1 and 2. In order to point out this "degeneration", a new parameter α is introduced, leading to the FEL G(n, α, f ). At a given n, we must calculate the (restricted) partition sum over the configurations α to get G(n, f ). In order to find G(n, α, f ), we consider the following

G(n, α, f ) = G form (n, α) + G ss str (n, α, f ) -f ∆x ss l , (8)
where G form (n, α) is the free energy of formation of the configuration (n, α), G ss str (n, α, f ) is the strain energy of the unfolded ssDNA and f ∆x ss l is a Legendre term (for the definition of ∆x ss l see eq. ( 10) below). The free energy of formation G form (n, α) is written as

G form (n, α) = all bps g µ,µ+1 + G junc (n) + G (1) loop (n) + G (2) loop (n) . (9)
The first term refers to the nearest neighbour model developed in [33,34], which specifies the interaction g µ,µ+1 between a base pair µ and the directly adjacent one µ+1. It was shown to provide reasonable agreement with experiments [24,29,35]. For example, applying this model onto a sequence 5 -TCCAG. . . -3 and its complementary part 3 -AGGTC. . . -5 , the stack energy reads G stack = g TC/AG + g CC/GG + g CA/GT + g AG/TC + . . . The most recent values of g µ,µ+1 lie in the range of -2.37 to -0.84 kcal/mol at 25°C [27]. The terms G junc (n), G (1) loop (n) and G (2) loop (n) in eq. ( 9) describe the free energy reduction due to the release of the stem-hairpin-junction and end-loops and are estimated from [START_REF] Santalucia | Proc. Natl. Acad. Sci[END_REF][START_REF] Zuker | [END_REF] as G (1) loop = G (2) loop = 1.58 kcal/mol and G junc = 4.90 kcal/mol. The strain energy G ss str (n, α, f ) of the unfolded single-stranded part [24] with contour length l = l(n, α)1 can be calculated from the work needed to stretch the unpaired bases. We like to remind the reader that we denote the elongation of the mean end-to-end distance of the ssDNA in force direction by u l ( f ). In what follows we chose to work with the FJC model, see eq. ( 4), with parameters of [27] instead of using the WLC model, despite the fact that both approaches give similar good results (cp. table 1). Since u l ( f ) is monotonously increasing with f , it has an inverse f l (u) = u -1 l ( f ), which is the force that is exerted by a ssDNA chain with contour length l, if its mean end-to-end distance is elongated by u. Accordingly, setting ∆x ss l = u l ( f ), we can write

G ss str (n, α, f ) = ∆x ss l 0 du f l (u ) = f ∆x ss l - f 0 d f u l ( f ) . (10)
Finally, we computed G(n, f ) by

G(n, f ) = -k B T ln α exp - G(n, α, f ) k B T . (11)
In fig. 7(b) the FEL is depicted for f = 0 and f = 16.21 pN. It exhibits the behaviour anticipated in fig. 7(a): for zero force it has minima at the stable/metastable states and it becomes tilted with rising force. When approaching the force regime where the rips occur in fig. 4(a), the levels of the minima become comparable. We want to point out that the DNA sequences shown in fig. 1 have been designed deliberately to yield the multiplestate structure seen in fig. 7(b). This gives us some confidence in the model underlying the construction of the G(n, α, f ) in eq. (8). Based on the FEL we can easily calculate the transition probability W( f i | f i-1 ) for the first rupture force f i in the ith transition if the first rupture force was f i-1 in the (i -1)th transition. The result is

W( f i | f i-1 ) = γ 0 i r γ i,i+1 ( f ) exp - γ 0 i r f f * i-1 d f γ i,i+1 ( f ) , (12)
where f * i = f i -∆ f i (and f 0 = f * 0 = 0); r was the loading rate, see sec. 2, and γ 0 i and γ i,i+1 ( f ) were defined in eq. ( 6). The activation energy ∆G i,i+1 ( f ) appearing in eq. ( 7) was calculated, as indicated in fig. 7(a), from the G(n, f ) by determining the energy G min i ( f ) of the local minimum belonging to state i and the saddle point energy G saddle i,i+1 ( f ) of the ith transition between the ith and (i+1)th state, ∆G i,i+1 ( f ) = G saddle i,i+1 ( f )-G min i ( f ). The attempt rate γ 0 i was used as the only fitting parameter. For the joint probability density of the three first rupture forces we then obtain

Ψ 3 ( f 1 , f 2 , f 3 ) = W( f 1 |0) W( f 2 | f 1 ) W( f 3 | f 2 ) , (13)
which allows us to calculate the distributions shown in fig. 6(a). Figure 8 displays the histograms for the three rips from fig. 6(a) in comparison with the distributions calculated from our theory. In view of the available statistics (82 cycles, see sec. 2), the agreement is quite satisfactory.

Conclusions

An important class of biophysical studies is the investigation of junctions in molecules since they present manifold ways to interact with other substances, for instance cations. Three-way junctions are especially interesting because metal ions such as magnesium [36,37] can bind to them and alter the tertiary structure. Here a first step of such a study is presented where we investigate a molecule with a three-way junction alone, without in vivo relevant substances.

One of our aims of this work was to study whether the construction and kinetics of more complex DNA molecules with richer folding-unfolding behaviour can be described by proper extensions of theories developed successfully for two-state systems so far. Our results show that this is indeed possible, at least for the unfolding trajectories. A triplebranch molecule has been specifically designed to produce a four-state system based on a model for the free energy landscape. This design was successful and we were able to prove the existence of these states by the emergence of associated force rips in pulling experiments. The first rupture forces have been systematically recorded in these pulling experiments and their distributions have been calculated. A transition rate theory based on the free energy landscape was successful in describing these distributions.

Two patterns have been found in the measured unfolding trajectories, one indicating the anticipated unfolding behaviour and the other one pointing to the occurrence of irreversible molecular fraying. This characterisation was possible by connecting the extracted force jump values to the change in the number of opened bps at each transition. For this estimation we tested the validity of two polymer models for the elastic response of ssDNA (FJC and WLC) and different sets of parameters in order to find the best agreement with the expected values. This analysis is useful to compare the elastic properties measured in DNA unzipping experiments with those obtained by stretching ssDNA polymers [27].

One class of molecules required a smaller force than anticipated to unfold the stem since some bps at its basis are partly or completely melted due to photodamaging. Permanent molecular fraying is an usually undesired, but frequent effect in single-molecule studies and deserves special attention in order to reduce its distorting influence on experimental results. It is important to find means to avoid irreversible fraying since not fully closed molecules change the measured unfolding-folding trajectories so that an average over all molecules, including permanently frayed ones, can lead to deviations of the real values and to misinterpreted results. To improve the statistics of the results, it is therefore necessary to identify irreversibly frayed molecules and remove them from the analysis. Within our analysis, we found a useful method to identify permanent molecular fraying. With the appropriate parameters for the polymer models one can estimate the number of irreversibly frayed bps at the basis of the stem. Further experiments performed with silica beads instead of polystyrene microspheres could clarify under which conditions permanent molecular fraying can be decreased. All this knowledge paves the way to further interesting studies such as the refolding problem [28] or two topics which require further experimental research. Firstly, the binding of metal ions to the three-way junctions could be examined in order to find out about structural changes due to the formation of tertiary contacts and therefore altered kinetics of the unfolding process [37]. Secondly, the translocation motion of helicases that unwind dsDNA could be addressed as another interesting subject, including the investigation of how they move along bifurcation points. Eventually, the kinetic approach for the prediction of the probability distributions of the first rupture forces could be a well suited starting point for future in-depth modelling in this domain on a more microscopic basis. folding unfolding

F

d ′ 0 x(n, f ) = u l (f ) + d ′ 0 n = 21 x(n, f ) = d 0 d 0 x(n, f ) = u l (f ) x(n, f ) = u l (f ) + d 0 d 0 d 0 n = 37 n = 53 l = 2nd l = 2nd + n loop d l = 2nd + 2n loop d n = 0
f 1 f 3 f 2 f 1,rf f 2,rf f 3,rf (b)

4 .

 4 Different unfolding patterns related to the number of opened bps

  )Here Y denotes the Young modulus, b is the Kuhn length, k B the Boltzmann constant and T the temperature. Typical values of the Kuhn length and the Young modulus under working conditions of T 25 °C and 1 M NaCl concentration are b = 1.42 nm and Y = 812 pN[18] or, as published recently, b = 1.15 nm and Y = ∞[27]

Figure 3 .

 3 Figure 3. The four stable or metastable states of the triple-branch molecule: 1 -folded molecule, 2 -unfolded stem, 3stem and one hairpin branch unfolded, 4 -unfolded molecule. The values of x(n, f ) refer to the end-to-end distance given in eq. (2) with the number n of opened bps corresponding to the molecular construct shown in fig.1and the contour length l according to eq. (3).

Figure 4 .

 4 Figure 4. Force as a function of change in the trap-pipette distance (a) for one typical unfolding trajectory for the first investigated molecule and (b) for one typical unfolding and folding trajectory for the second investigated triple-branch molecule. Indicated are the three first rupture forces belonging to the transitions 1 ( f 1 ), 2 ( f 2 ) and 3 ( f 3 ) and the respective refolding forces, labelled f i,rf . The inset in (a) shows four further unfolding curves for the first investigated molecule.

Figure 5 .Figure 6 .

 56 Figure 5. The first rupture forces and force jump values of the unfolding trajectories shown in fig. 4 are extracted as indicated here from the intersection of the fitted grey lines and the transitions (vertical lines). In part (a) the procedure is shown for the molecule used in the theoretical analysis in sec. 5 and in (b) for the molecule with permanent fraying behaviour.

Figure 7 .Figure 8 . 1 = 9 • 10 6 Hz, γ 0 2 = 1 . 1 • 10 6

 7819621110 Figure 7. (a) Sketch of the expected FEL for the triple-branch molecule at zero force as a function of the number n of opened bps. (b) FEL calculated from eqs. (8) to (11) at two different forces in units of the thermal energy k B T .

  and 21 bps n loop d, n =21 + 16 bps 2n loop d, n =21 + 16 + 16 bps

	,

Table 1 .

 1 Overview over the change in the number of opened bps for different models and parameters for the two representative molecules. The numbers in brackets are the standard deviations.

	and for the WLC model

  Figure2. Sketch of the experimental setup. In our experiments we measure the relative distance X rather than its absolute value.

	December 1, 2010 19:33 Philosophical Magazine engel˙etal˙complex˙systems˙12 December 1, 2010 19:33 Philosophical Magazine engel˙etal˙complex˙systems˙12 December 1, 2010 19:33 Philosophical Magazine engel˙etal˙complex˙systems˙12
	1	DNA triple-branch molecule 2 3	4
		hairpin branch 1 hairpin branch 2 Hairpin branches Loop Loop
		16bps + 4bs loop	16bps + 4bs loop
	Optical trap	Stem
				Micropipette
	AD bead	Handle 1		Handle 2	bead SA
	o r F o r F o r molecular extension absolute distance F o r
		P P P P
		e e r e e r e e stem 21bps e r e r
		R e R e R R e e	Total: 114bs
			v i e v i e v i e v i e
				w w w w
				O n l O n l O n l O n l
					y y y y

Figure 1. Structure of the triple-branch DNA molecule.

In our experiments we measure the relative distance between trap and pipette, X, rather than the absolute value. The force f exerted on the molecular construct leads to a displacement f /k b of the bead in the optical trap lead, where k b = 0.08 pN/nm is the rigidity of the trap. Hence, the relative distance X is related to the relative molecular extension x m (see fig.

2) byx m = Xf /k b .

Due to this simplification, the ∆n 1 of the first rip is likely to be slightly underestimated and the second rip's ∆n

overestimated. However, it will not affect the change in the total number ∆n tot of opened bps since we consider the change of x(n, f ) and the d 0 contributions will cancel each other out.

We like to note that the checking of the change in the number of opened bps can be, in principle, also applied to nonpermanent, reversible molecular fraying.

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

In previous works of some of the authors this contour length was denoted by l n,α to emphasize the dependence on n (and, in addition, α here). This dependence is caused by the change of the contour length in the transitions. For easier reading we suppress to give it explicitely in the following. Further details about the contour length were already discussed in sec. 4.

Acknowledgements

S. E. thanks the Deutscher Akademischer Austauschdienst (DAAD) for providing financial support (FREE MOVER and PROMOS) for stays at the Small Biosystems Lab in Barcelona where the experiments have been performed. A. A. is supported by grant AP2007-00995. F. R. is supported by the grants FIS2007-3454, Icrea Academia 2008 and HFSP (RGP55-2008).