Marco Ripà
email: marcokrt1984@yahoo.it

On some open problems concerning perfect powers

Keywords: Open problem, Perfect power, Perfect square, Conjecture, Integer sequence. 2020 Mathematics Subject Classification: 00A08 (Primary) 11B50, 11-04 (Secondary)

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In late 2010, the author of this paper found a recreative open problem by Kenichiro Kashihara (see [START_REF] Kashihara | Comments and topics on Smarandache notions and problems[END_REF], open problem number 30, p. 25) concerning the sequence A001292 of the On-Line Encyclopedia of Integer Sequences (OEIS) [START_REF] Muller | Sequence A001292 in the On-line Encyclopedia of Integer Sequences[END_REF]. Kashihara's problem number 30 consists of two independent parts and the author solved the first one quite easily at the time (the complete solution can be found in [START_REF] Ripà | On prime factors in old and new sequences of integers[END_REF], Section 3.3, pp. 12-15), since it asks to find the probability 0 < p(c) < 1 that the trailing digit of the generic term of the sequence A001292 is c ∈ 0, 1, 2, . . . , 9 and the formula provided in [START_REF] Ripà | On prime factors in old and new sequences of integers[END_REF] shows that p(c) = 11-c 55 for any c ̸ = 0, whereas p(0) = 0.018 (e.g., if c = 7, then p(7) = 4 55 = 0.072). In the present paper, we will focus ourselves on the second part of the above mentioned Kashihara's problem number 30, asking how many elements of the sequence A001292 are perfect powers, since Kashihara conjectured that there are none. Now, bearing in mind that a perfect power of an integer d > 0 is a natural number k ≥ 2 such that a k = d, where also a is a positive integer, we could point out that A001292(1) = 1 can be considered as a solution and argue how this disproves the conjecture, but (from here on) we will disregard this special case and assume that we are looking for a nontrivial counterexample to Kashihara's conjecture.

Lastly, Section 3 is to devoted to introducing a (quite improbable) conjecture concerning perfect powers belonging to the OEIS sequence A352991 [START_REF] Ripà | Sequence A352991 in the On-line Encyclopedia of Integer Sequences[END_REF][START_REF] Ripà | Sequence A353025 in the On-line Encyclopedia of Integer Sequences[END_REF].

The exclusion criterion

In order to be clear on the invoked OEIS sequences, let us introduce a few useful definitions.

Definition 1.

We define the m-th term of the OEIS sequence A007908 as A007908(m) := 123 . . . (m -1) m, where m ∈ Z + . Definition 2. We define the sequence A001292 of the OEIS as the concatenations (sorted in ascending order) of every cyclic permutation of the elements of the sequence A007908 (e.g., given m = 3, A001292(A007908(3)) = 123, 231, 312). Definition 3. We define the OEIS sequence A352991 as the concatenation of all the distinct permutations of the first strictly positive m integers, sorted into ascending order (e.g., 12345671089 is a term of the sequence, while 12345670189 does not belong to A352991, even if all the digits of the string 123 . . . 910 appear once and only once, since "10" is missed).

After having checked the first 100128 terms of the sequence A001292 (see Appendix), exploring any exponent at or above two, we have not found any perfect power so that Kashihara's conjecture has been verified up to 10 1235 (i.e., the 100129-th term of A001292 is the smallest cyclic permutation of A007908(448) and is greater than 10 1235 by construction). Moreover, we can prove the following Theorem 1, concerning the sequence A352991 which includes every term of A001292.

Theorem 1. For any m > 1, A352991(n) cannot be a perfect power of an integer if A352991(n) is a permutation of A007908(m) and m : m ≡ {2, 3, 5, 6}(mod 9).

Proof. By definition, A007908(m) [START_REF] Muller | Sequence A007908 in the On-line Encyclopedia of Integer Sequences[END_REF] cannot be a perfect power if 123 . . . (m -1) m is divisible by 3 and it is not divisible by 3 2 . Thus, from the well-known divisibility by 3 and 9 criteria, m :

(3 | m j=1 j) ∧ (3 2 ∤ m j=1 j
) is a sufficient, but not necessary, condition for letting us disregard any permutation of 123 . . . (m -1) m (i.e., given m, if a generic permutation of A007908(m) is divisible by 3 and is not congruent to 0(mod 9), then all the permutations of A007908(m) are divisible by 3 once and only once, since the commutativity property holds for addition).

It follows that, for any n ∈ Z + , A352991(n) cannot be a perfect power if it is a permutation of the string 123 . . . (m -1) m, where m is such that A134804(m) is divisible by 3. Therefore, the residue modulo 9 of every perfect power belonging to A352991 cannot be 2 or 3 or 5 or 6, and this concludes the proof of Theorem 1.

Corollary 1. Kashihara's conjecture is true for the concatenation of any cyclic permutation of A007908(m), where m : m ≡ {2, 3, 5, 6}(mod 9) ∨ m < 448.

Proof. We observe that A001292 is a subsequence of A352991 [START_REF] Muller | Sequence A001292 in the On-line Encyclopedia of Integer Sequences[END_REF][START_REF] Ripà | Sequence A352991 in the On-line Encyclopedia of Integer Sequences[END_REF]. By invoking Theorem 1, we can state that every perfect power candidate has to be the concatenation of a (cyclic) permutation of A007908(m), where m is such that m : m ≡ {2, 3, 5, 6}(mod 9). On the other hand, all the remaining terms up to 99 100 101 . . . 445 446 447 1 2 3 . . . 96 97 98 have been directly checked (see Appendix for details) and no perfect power has been found.

Therefore, Corollary 1 confirms Kashihara's conjecture for any term of A001292 such that m is congruent to {2, 3, 5, 6}(mod 9) or m ≤ 447.

Corollary 2. ∄n : A353025(n) ≡ {2, 3, 4, 5, 6, 7, 8}(mod 9), and any term of A001292 cannot be a perfect power if its digital root is not equal to 0 or 1.

Proof. Trivially, 10 ≡ 1(mod 9) and also (1 + 0) ≡ 1(mod 9) so that any positive integer is congruent modulo 9 to its digital root. Now, we observe that any terms of A001292 belongs to A353025.

Since from Theorem 1 it follows that every term of the sequence A353025 [START_REF] Ripà | Sequence A353025 in the On-line Encyclopedia of Integer Sequences[END_REF] is a special permutation of A007908(m) which is characterized by m ≡ {0, 1, 4, 7, 8}(mod 9), in order to prove Corollary 2, it is sufficient to note that

m j=1 j ≡                    0(mod 9) if m : m ≡ 0(mod 9) 1(mod 9) if m : m ≡ 1(mod 9) 1(mod 9) if m : m ≡ 4(mod 9) 1(mod 9) if m : m ≡ 7(mod 9) 0(mod 9) if m : m ≡ 8(mod 9)
.

(

) 1
Remark 1. A well-known property of integers is that every perfect power which is congruent modulo 5 to 0 is also necessarily congruent to {0, 25, 75}(mod 100), while if a perfect power is congruent modulo 10 to 6, then its second last digit is odd. Thus, we are free to combine these additional constraints with Corollary 2 in order to reduce the number of perfect power candidates among the terms of A352991.

Perfect cubes in A353025

In the first half of April 2022, playing with Kashihara's conjecture, a more risky (very likely false but really hard to disprove by brute force) conjecture arose, it is as follows.

Conjecture 1. Let n ∈ N -{0, 1} be given. We (provocatively) conjecture that if n is such that A352991(n) is a perfect power of an integer, then ∄k ∈ N -{0, 1, 2} :

A352991(n) = c k , c ∈ N.
On April 16, 2022, a direct search was performed by the author on the first 10 7 terms of the sequence and no counterexample has been found (42 perfect squares only). A few days later, Aldo Roberto Pessolano, performed a deeper search running the Mathematica codes published in Appendix, without finding any counterexample and thus confirming Conjecture 1 (at least) up to the smallest permutation of A007908(22) (i.e., for any term of A352991 which is strictly greater than 1 and smaller than 10111121314151617181920212223456789).

Remark 2. If confirmed, Conjecture 1 would imply that all the perfect powers (stricly greater than 1) in A352991 are perfect squares and nothing more (no cube, no square of square, and so forth). Nevertheless, under the (arbitrary, but perfectly reasonable) assumption of a standard probability distribution of the cubes in A352991 (i.e., we are assuming that |n∈N:A352991(n)≤m∧A352991(n)

1 3 | |n∈N:A352991(n)≤m| ∼ = |n∈N:A000578(n)≤m| m
holds for any sufficiently large m ∈ N), we would guess the existence of infinitely many counterexamples to Conjecture 1, even if the smallest one is expected to occur in the interval [10 58 , 10 65 ]. On the other hand, the same argument would corroborate Kashihara's conjecture, since the number of perfect powers belonging to A001292 cannot probabilistically exceed

2 • +∞ k=2 +∞ j=309 4 • j • (10 4•j k -10 4•j-1 k ) 9 • 10 (4•j-1) ≈ 8 9 • +∞ j=309 j • (10 4•j 2 -10 4•j-1 2 ) 10 (4•j-1) ≈ 0. (2)
Additional open problems. Does the sequence A353025 have infinitely many perfect squares, infinitely many perfect cubes, infinitely many perfect squares of squares, and so forth? Which is the smallest nontrivial perfect cube (if any) belonging to A353025 (we point out that all the terms greater than one and below (1.01 • 10) 40 have been checked without finding any cube)?

Conclusion

Kashihara's open problem number 30 has not been completely solved yet. Even if the first part, concerning the probability that the trailing digit of A001292(n) is c = 1, 2, . . . , 9, was solved by the author a dozen of years ago [START_REF] Ripà | On prime factors in old and new sequences of integers[END_REF], the second part and the related conjecture still needs a proof or a nontrivial counterexample (the smallest candidate has 1236 digits). Moreover, in the present paper, we have introduced a wider speculation that allow us to ask to ourselves how to find a term of the OEIS sequence A353025 (disregarding A353025(1)) which is not a perfect square; a challenging open problem, considering that there is not any perfect cube among the terms of A352991 in the interval (1, 10 40 ]. On the other hand, the following code returns the complete list of the smallest 42 perfect squares belonging to A352991.

z = 1; Do[r = Range[k]; n = ToExpression[StringJoin[ToString[#]&/@r]]; If[And[Mod[n, 9] != 3, Mod[n, 9] != 6], d = DigitCount[n]; s = IntegerPart[Sqrt[10ˆ(IntegerLength[n] -1)]]; f = IntegerPart[Sqrt[10ˆ(IntegerLength[n])]]; Do[y = xˆ2; If[DigitCount[y] == d, c = True; Do[ If[Not[StringContainsQ[ToString[y], ToString[i]]], c = False ], {i, 10, k}]; If[c, Print[z, " ", y]; z++] ], {x, s, f}] ], {k, 2, 10}]
These 42 perfect squares correspond to all the perfect powers in (1, 10 16 ] belonging to A352991, while the next perfect square is 10135681742311129 (we observe that 100676123 2 is a permutation of 123 . . . 16, as suggested by the statement of Theorem 1). In the end, our tests have finally confirmed that all the perfect powers which are smaller than 10 34 and that belong to the OEIS sequence A352991 are perfect squares (only). At the present time, Conjecture 1 has been tested for every integer smaller than 10111121314151617181920212223456789, and no counterexample has been found yet.

  Aldo Roberto Pessolano helped the author of the present paper by verifying Kashihara's conjecture and Conjecture 1 for a very large number of terms. All the provided Mathematica codes run on the M1 processor of his Apple MacBook Air (2020). Kashihara's conjecture has

	z = True;
	Break[]]&/@{2, 3, 5, 6, 7};
	h++;g = 2ˆh]];
	If[z, Break[],
	Print["1..", k, " checked."]],
	{k, 2, 21}]
	5 Appendix

Acknowledgments

We sincerely thank Aldo Roberto Pessolano for having helped us very much with the search of perfect powers belonging to A001292 and A352991, letting us confirm Kashihara's conjecture up to 10 1235 and Conjecture 1 up to (1.01 • 10) 34 .

been currently tested up to the 100128-th term of A001292 and we confirm that it holds for every element of the set A001292(2), A001292(3), . . . , A001292(100128).

The search reached the term 99 100 . . . 446 447 1 2 . . . 97 98 ≈ 9.91 • 10 1232 in 28823 seconds (about 8 hours of calculations) and the code is as follows:

About our investigation on the perfect powers in A352991, Pessolano has recently completed the direct check of every term of A352991 which falls in the interval (1, 987654322120191817161514131211110] (see the code below). As expected, the test has not returned any perfect power above two.

, #, d, k], Print[x, "ˆ", #, " = ", y]; z = True; Break[]],{x, s, f}])&/@p;g = 2ˆh; While[g < n,If[q[#, h, d, k], Print[x, "ˆ", h, " = ", y];


Eric Guiot
email: guiot.eric_1@yahoo.fr

Trajectories of two-dimensional harmonic oscillators in rotating frame: application to Foucault pendulum problem and centered trochoid trajectories

Keywords: Centered trochoid curves, harmonic oscillator, non inertial frame, Simple gravity pendulum

The motion of a classical point particle submitted to a central restoring force in a noninertial frame is studied. Exact trajectories are provided and several classical cases are detailed. Results are applied to the classical case of the Foucault pendulum problem, and centered trochoid trajectories are highlighted. These original solutions are linked to concerns of contemporary physics. .

Introduction

It is well known that studies of the trajectories of particles in non-inertial frames are major concerns of physics. Among them, behaviors of harmonic oscillators in such frames are particularly interesting, given the many possible contemporary physical applications. We are thinking, for example, of studies concerning spectroscopy or plasma physics, whether in a classical or quantum framework. Another major application, which will be particularly studied in this article, concerns the Foucault pendulum problem. Indeed, this problem is generally studied in the small-angle approximation, which implies that the pendulum becomes, in this case, a two-dimensional harmonic oscillator.

The purpose of this paper is to contribute to this topic, by presenting an original study of the trajectory of these harmonic oscillators in non-inertial frames, more precisely in time varying rotating frames. Indeed, using the laws of classical mechanics, we will present original solutions, which seem to be known only for certain limit cases. The particular case where the speed of rotation of the non-inertial frame is constant will be highlighted, and we will show that the corresponding trajectories are the centered trochoid curves. In a second part, we will apply these general results to the case of Foucault's pendulum. This will allow us to find certain limiting cases already known, but also to characterize geometrically several possible trajectories and to determine the laws of motion.

There are many reasons for us to apply our results to this particular problem. First of all, it is, of course, one of the most important popular experiment in the history of Physics, which was presented for the first time by the French physicist Foucault to the Parisian public in 1851 [1], and was subsequently reproduced in many museums around the world. Schematically, the experiment consists in setting in motion a long and heavy pendulum and measuring the slow rotation of the plane of oscillation, due to the Earth rotation. It is then possible to highlight this rotation and to deduce, for example, the latitude of the pendulum. Naturally, studies of the behavior of the pendulum has been presented in many physics text books and various solutions have been presented (the reader can consult, for example, the classical references [START_REF] Taylor | Classical Mechanics[END_REF][START_REF] Gregory | Classical Mechanics: An Undergraduate Text[END_REF][START_REF] Lamb | higher mechanics[END_REF]). Unfortunately, these works, and more particularly the establishment of equations of motion and determination of geometrics characteristics of trajectories, do not seem to be unanimous, and therefore remains an important concern of the contemporary Physics. This explains why research papers are still regularly published on this topic with the aim, in general, to provide the most precise possible laws of motion.

Another reason concerns the possible practical applications of this experiment. For example, it is planned for measurement of the lense-Thirring effect and, consequently to improve our knowledge of the theories of Gravitation and of the space time. Moreover, for the reason the extension of the pendulum principle allows the development of gyroscopic systems, very numerous applications in many fields of applied physics and technology are considered.

Therefore, paper is organized as follows: in a first part, we study the behavior of a two-dimensional harmonic oscillator in all kind of time varying rotating frame. Secondly, we apply our results to the particular case of the Foucault pendulum.

The reader will therefore find in this article progress concerning several areas of the physics: on the one hand, from a geometric point of view, a new approach to trochoid curves. On the other hand, exact solutions of the movement two dimensional harmonic oscillator in time varying rotating frame and application to the Foucault's pendulum problem, in the approximation of small angles. To conclude, we attract his attention that this work has been carried out within the framework of classical mechanics (quantum and relativistic effects have been neglected).

Two dimensional harmonic oscillator in time varying rotating frame 2.1 Present the problem

(Note, to simplify the calculation, the masses were, in this part, considered equal to 1, so that the forces and the accelerations are confused)

We begin the study by introducing a Galilean system of coordinate (𝑂, 𝑋 0 ⃗⃗⃗⃗ , 𝑌 0 ⃗⃗⃗ , 𝑍 0 ⃗⃗⃗⃗ ) where a point particle 𝑀 moves under the influence of a restoring force centered at 𝑂 and given with

𝐹 = -𝑘𝑂𝑀 ⃗⃗⃗⃗⃗⃗ = -𝑘𝑟 = -𝑤 2 𝑟 ( 1
)
Where 𝑤 is the pulsation of the pendulum. We consider now a second system of coordinate (𝑂, 𝑥 , 𝑦 , 𝑧 ). This frame is animated by a rotation movement around the z axis such

Ω ⃗⃗ (𝑡) = Ω𝑧
And is consequently non-inertial. As we presented it, goal of the paper is to describes the trajectory of 𝑀 in this frame.

To obtain the equation of motion we have to consider the fictitious forces, defined, according to the laws of classical mechanics, by

𝐹 𝐹𝑖 ⃗⃗⃗⃗⃗ = 𝑎 0 ⃗⃗⃗⃗ -Ω ⃗⃗ × Ω ⃗⃗ × 𝑟 -2Ω ⃗⃗ × 𝑑r 𝑑𝑡 - 𝑑Ω ⃗⃗ 𝑑𝑡 × 𝑟 (2)
Where we successively recognize the acceleration of translation between the two frames, centrifugal force, the force of Coriolis and of Euler (see for example classical textbooks in references [START_REF] Taylor | Classical Mechanics[END_REF][START_REF] Gregory | Classical Mechanics: An Undergraduate Text[END_REF]).

Respecting our conditions, the translational acceleration is zero and the equation of the dynamics becomes in this frame

𝐹 = -𝑘𝑟 -Ω ⃗⃗ × Ω ⃗⃗ × 𝑟 -2Ω ⃗⃗ × 𝑑r 𝑑𝑡 - 𝑑Ω ⃗⃗ 𝑑𝑡 × 𝑟 (3)

Solving

Equation of motion leads to the system :

{ 𝑋 ̈= -𝑤 2 𝑋+Ω 2 𝑋 + 2Ω𝑌 ̇+ 𝑌Ω Ẏ̈= -𝑤 2 𝑌+Ω 2 𝑌 -2Ω𝑋 ̇-𝑋Ω ̇} (4)
Again, these equations are not original (see for example [START_REF] Taylor | Classical Mechanics[END_REF]). The classical method of solving consists in introducing a complex number such that

𝑢 = 𝑋 + 𝑗𝑌 Differential equation (4) becomes 𝑑 2 𝑢 𝑑𝑡 2 = -(𝑤 2 -Ω 2 )𝑢 -2𝑗Ω 𝑑𝑢 𝑑𝑡 -𝑗Ω ̇𝑢 (5)
To solve, we choose to investigate a solution in the original following form:

𝑢 = 𝑟 1 𝑒 𝑗(𝑤𝑡+𝑓(𝑡))
And (5) becomes

𝑗𝑓 ̈-(𝑤 + 𝑓 ̇)2 = -𝑤 2 + Ω 2 + 2Ω(𝑤 + 𝑓 ̇) -𝑗Ω ̇
Identifying the imaginary and real parts we get the system of equations :

{ 𝑓 ̈= -Ω ̇ 6.1 -(𝑤 + 𝑓 ̇)2 = -𝑤 2 + Ω 2 + 2Ω(𝑤 + 𝑓 ̇) 6.2 }
By considering the constants of integration as null (6.1) leads to

𝑓(𝑡) = -∫ Ω(t)𝑑𝑡 (7)
Introducing 𝑓 ̇= -Ω In 6..2 we verify easily the solution. Moreover, we note that a second solution exits, given by

𝑢 = 𝑟 2 𝑒 𝑗(-𝑤𝑡+𝑓(𝑡))
Because it leads to the same relation ( 7) and also satisfies real relation, given this time with

-(-𝑤 + 𝑓 ̇)2 = -𝑤 2 + Ω 2 + 2Ω(-𝑤 + 𝑓 ̇)
Whereas these results we obtain thus a general solution such

𝑢 = 𝑟 1 𝑒 𝑗(𝑤𝑡+𝑓(𝑡)) + 𝑟 2 𝑒 𝑗(-𝑤𝑡+𝑓(𝑡))
Rewritting the expression, we obtain :

{ 𝑋 = (𝑟 1 + 𝑟 2 ) cos 𝑤𝑡 cos f(t) -(𝑟 1 -𝑟 2 )sin 𝑤𝑡 sin f(t) 𝑌 = (𝑟 1 + 𝑟 2 )cos 𝑤𝑡 sin f(t) + (𝑟 1 -𝑟 2 )sin 𝑤𝑡 cos f(t) } (8)
Now, consider the limiting cases for which we know the solution. This is, of course, the case where Ω = 0. Then, the point particle describes naturally an ellipse of Lissajous-Bowditch, define (with correct initial conditions) by

{ 𝑋 = 𝑎 cos 𝑤 𝑡 𝑌 = 𝑏 sin 𝑤 𝑡 }
Where (𝑎, 𝑏) are the semi major and minor axis of the conic. We deduce that

{ 𝑟 1 + 𝑟 2 = 𝑎 𝑟 1 -𝑟 2 = 𝑏 }
Finally the solution can be written in the form of simple trigonometric functions, for example

{ 𝑋 = 𝑎 cos 𝑤𝑡 cos f(t) -𝑏 sin 𝑤𝑡 sin f(t) 𝑌 = 𝑎 cos 𝑤𝑡 sin f(t) + 𝑏 sin 𝑤𝑡 cos f(t) } (8)
Note that it is possible to verify these results, by calculating the successive derivatives of these solutions and introducing them into the equation of motion [START_REF] Lamb | higher mechanics[END_REF]. Also note that these solutions seems original in general cases.

Presentation of trajectories

E.Guiot

Trajectories of two-dimensional harmonic oscillators in rotating frame: application to Foucault pendulum problem and centered trochoid trajectories Note that these trajectories are not unknown but represent a part of the centered trochoid curves (more exactly the hypotrochoid curves) as defined, for example, in reference [5].

Among these trajectories, several particular cases can be noticed : an emblematic case concern the rose trajectory, corresponding on a mono dimensional oscillator (𝑏 = 0). In this case, the angular speed is constant and simply equal to Ω [START_REF] Guiot | Two-dimensional Harmonic oscillator and centered Trochoïde curves: Trajectories under constant magnetic field[END_REF]. Others well known trajectories can be reobtain, such the hypocycloid, etc.. Certains of them will be drawn in the following parts.

General case

It is also possible to study the trajectories of the harmonic oscillator in a time-dependent rotating frame. Consider for example the classic example Ω ⃗⃗ = Ω(1 + 𝛽 cos 𝜈𝑡)𝑧 We obtain using (7) the expression

𝑓(𝑡) = -Ω(𝑡 + 𝛽 𝜈 sin 𝜈𝑡)
In figures 2.1 and 2.2 we present the modifications of the trajectories due to the addition of a timedependent part to the constant rotational speed Ω Consider also the case of the monodimensional oscillator (figure 2.3) Another classic case concerns a speed of rotation which gradually decreases until it cancels out: Considering for example the case Ω = Ω 0 𝑒 -𝜆𝑡 . We obtain the trajectories this time become the trochoids corresponding to the constant rotation speed.

Application to Foucault pendulum problem

As we introduced, we now apply these theoretical results to the case of the Foucault pendulum problem. First, we quickly present the classical case of the simple pendulum in an inertial frame.

Simple gravity pendulum

Simple pendulum is naturally a well-known oscillating system, which has often be described in classical text books (see for example [START_REF] Taylor | Classical Mechanics[END_REF]). In figure 3, the force diagram and the system of coordinate is detailed. Note that friction forces and mass of the rod are neglected, and that length 𝑃𝑀 = 𝐿 is considered constant. The system of coordinate (𝑂, 𝑥 , 𝑦 , 𝑧 ) is, in this case, Galilean (because Earth 's rotation is neglected). We can write the second equation of dynamics, whic leads to the following relation

𝑃 ⃗ + 𝑇 ⃗ = 𝑚𝑎 ( 10
)
Where 𝑎 is the acceleration of the point mass 𝑀(𝑚). In the small angle approximation, the motion is considered planar and contained in the plane (𝑂, 𝑥 , 𝑦 ). (10) become thus

{ -𝑇 sin 𝜃 = 𝑚𝑦T cos 𝜃 -𝑚𝑔 𝑂 = 𝑚𝑧̈} (11)
-The centrifugal force (see reference [START_REF] Taylor | Classical Mechanics[END_REF]) is directed toward the rotation axis of the Earth and is independent of the motion of the particle in the rotating frame. Its major effect is to modify the gravitational acceleration, which becomes the effective gravitational acceleration given by

𝑔 = 𝑔 0 ⃗⃗⃗⃗ -𝑚Ω ⃗⃗ × (Ω ⃗⃗ × 𝑟 )
For this reason, the effect of this force on the pendulum's precession is usually neglected. It is the case in many theoretical Mechanics textbooks.

-The last force (Coriolis force) is thus considered as the only one which acts on the pendulum precession. In our system of coordinate, it can be written

𝐹 𝐶 ⃗⃗⃗⃗ = -2𝑚Ω ⃗⃗ × (𝑥ẋ + 𝑦̇𝑦 + 𝑧ż )
However, with the small angle approximation, we consider that 𝑧̇= 0 and we obtain thus the following expression in the plane (𝑂, 𝑥, 𝑦)

𝐹 𝐶 { 2𝑚𝑦̇Ω sin 𝜑 -2𝑚𝑥Ω sin 𝜑 }
Finally, the equations of motion can be written, in the rotating frame, by the system

{ 𝑥̈= - 𝑔 𝐿 𝑥 + 2Ω sin 𝜑 𝑥ẏ ̈= - 𝑔 𝐿 𝑥 -2Ω sin 𝜑 𝑦̇} (12)
Because of the results presented in the part 2.2, it makes sence to write the solutions in the form

{ 𝑥 = 𝑎 cos 𝑤𝑡 cos 𝑓𝑡 -𝑏 sin 𝑤𝑡 sin 𝑓𝑡 𝑦 = 𝑎 cos 𝑤𝑡 sin 𝑓𝑡 + 𝑏 sin 𝑤𝑡 cos 𝑓𝑡 } ( 13
)
Where 𝑓 is the pulsation given by 𝑓 = Ω sin 𝜑 And where the frequencies are linked with

𝑔 𝐿 = 𝑤 2 -𝑓 2
Trajectories are thus hypotrochoid curves. Depending on initial conditions, we can distinguish several classical cases.

Hypocycloid trajectories

E.Guiot

Trajectories of two-dimensional harmonic oscillators in rotating frame: application to Foucault pendulum problem and centered trochoid trajectories 12

To carry out the experiment, generally the pendulum is moved away from its equilibrium position on the x-axis, with a small angle 𝜃 0 and released at 𝑡 = 0 without initial velocity. These conditions can be written as follows :

{ 𝑎 = 𝐿𝜃 0 𝑥̇= 𝑦̇= 0 }
Introducing them in equation ( 13), we deduce the following relation :

𝑏 = -𝑎 𝑓 𝑤
To simplify the calculation we introduce the ratio

𝑟 1 𝑤 1 = 𝑟 2 |𝑤 2 |
We deduce that the trajectory is one of the hypocycloïd curves, among hem astroïd, deltoïd, etc.. see for example [7]. An example of these curves is provided in Figure 3.2 : Mechanical energy is given by

𝑝 = 2𝜋𝑛
𝐸 𝑀 = 1 2 𝑎 2 𝑤 2 (1 -𝑛 2 )

Rose trajectories

Consider the case where initial conditions are the following: at the origin of the time the pendulum is located at its equilibrium position, i.e. at the origin of the rotational frame. An impulsion is transmetted along the x-axis, such the initial velocity is 𝑉 0 ⃗⃗⃗ (0, 𝑦0, 0). As previously, we can obtain the expression of mechanical energy, given by

𝐸 𝑀 = 1 2 𝑏 2 𝑤 2
The interest of this trajectory is that it is the only one, among the wide variety of hypotrocoïd curves, which implies that the pendulum passes through the origin of the system of coordinate, at each oscillation. The time interval between each passage to O is given by the following relation

∆𝑇 = 𝜋 𝑤 Which leads to ∆𝑇 = 𝜋 √ 𝑔 𝐿 + 𝑓 2
Pulsation 𝑓 is thus obtained

𝑓 = √ [ 𝜋 ∆𝑇 ] 2 - 𝑔 𝐿
Therefore, it might be simpler and more accurate to measure this time interval that to measure the precession angle, as ikt usually seems to be.

General case

Discussion

Consider the general solutions presented in this article, about the behavior of harmonic oscillators in time varying rotating frames. First of all, it seems that these solutions are entirely original. Naturally, these trajectories can be applied to the Foucault pendulum only as a limiting case, if the experiment takes place classically, ie on Earth. Indeed, in this case the Euler force is zero. However, these trajectories can correspond to the movement of such a pendulum located on a time varying rotating frame, theoretically animated by all kinds of variable rotational speeds (we have presented some of them). For this reason, it is possible that these solutions are interesting to improve experiments, for example to verify gravitational effects, as we have introduced. Indeed, we know that the classic experiments would be very long and complex.

Moreover, these results can be applied to other areas, because they describe the theoretical motion of harmonic oscillators in rotating frames. Therefore, it is possible that these solutions can be applied, for example, in the study of gyroscopic movements, or such as those mentioned in the introduction to the article, or others.

Conclusion

We have detailed exact trajectories of two-dimensional harmonic oscillator in time varying rotating frame. In a second time, we have applied these results to the Foucault Pendulum problem. In particular, we suggest that trajectories are hypotrochoid curves, that seem a second original result. A simple experiment consists in supplying an impulse in the plane of the oscillations, and in verifying that the trajectory is well a rose curve. Later works could extend these results to other problems, such, for example, gyroscopic motion.

5 2 . 3 . 1

 231 Particular case: Two dimensionnal oscillator in constant speed of rotation : hypotrochoid curves Solutions becomes simply { 𝑥 = 𝑎 cos 𝑤𝑡 cos Ω𝑡 -𝑏 sin 𝑤𝑡 sin Ω𝑡 𝑦 = 𝑎 cos 𝑤𝑡 sin Ω𝑡 + 𝑏 sin 𝑤𝑡 cos Ω𝑡 } (9)

Figure 2 . 1 .Figure 2 . 2 .

 2122 Figure 2.1 . Left : hypocycloid of parameters : 𝑎 = 1, 𝑏 = -1/5, 𝑓 = 𝑡/5 Right : trajectory under parameter 𝑎 = 1, 𝑓 = 1/5𝑡 + 0.2 * sin(𝑡)

Figure 2 . 3 .

 23 Figure 2.3 . Left : rose curve of parameters : 𝑎 = 1, 𝑓 = 𝑡/5 Right : trajectory under parameter 𝑎 = 1, 𝑓 = 1/5(𝑡 + 0.2 * sin (𝑡)

Figure 2 . 4 .𝑓Figure 2 . 5 Figure 2 . 6

 242526 Figure 2.4 . Left : monodimensional oscillator : 𝑎 = 1, 𝑏 = 0, 𝑓 = -3 4 sin 𝑡

  Figure 3.1 Foucault Pendulum

  of the solutions are { 𝑋 = 𝑎[cos 𝑤𝑡 cos 𝑛𝑤𝑡 -𝑛 sin 𝑤𝑡 sin 𝑛𝑤𝑡] 𝑌 = 𝑎[cos 𝑤𝑡 sin 𝑛𝑤𝑡 + 𝑛 sin 𝑤𝑡 cos 𝑛𝑤𝑡] }

Figure 3 . 2

 32 Figure 3.2 𝑎 = 1, 𝑤 = 1, 𝑛 = 1/7

Figure 3 . 2

 32 Figure 3.2 Rose curve. 𝑎 = 1, 𝑤 = 1, 𝑛 = 7

Trajectories of two-dimensional harmonic oscillators in rotating frame: application to Foucault pendulum problem and centered trochoid trajectories

In this approximation we can also write following relations: Where 𝜃 0 is the initial angle and 𝑡 the time. Pulsation is given with

of this solution is that the Earth rotation is neglected and, consequently, that frame isn't exactly inertial: if we now consider this rotation, we have naturally to modify the equations of motion. This is the subject of the next chapter.

Physics: fictitious forces and trajectories

Thus consider this time the rotation of the Earth. The system of coordinate (𝑂, 𝑥 , 𝑦 , 𝑧 ) can no longer be considered Galilean but becomes, due to the Earth's rotation, a non-inertial rotating system frame. To detail corresponding fictitious forces we need to definite more precisely the system of coordinate (𝑂, 𝑥 , 𝑦 , 𝑧 ) that we use. Origin 0 is located at the latitude 𝜑 . To simplify the calculation, we choose the x-axis pointing east, y-axis north, and z-axis is vertical. 𝑣 (𝑥, 𝑦̇, 𝑧) is the speed of the pendulum and vector of rotation of the Earth is Ω ⃗⃗ . In this system, the coordinates of this vector is given by

We can now deduce the mathematical expression of the fictitious forces:

-In the classical problem, Euler's force is neglected because Earth's rotation speed is nearly constant.

E.Guiot

Trajectories of two-dimensional harmonic oscillators in rotating frame: application to Foucault pendulum problem and centered trochoid trajectories Consider this time that the pendulum is released such it would describe, if the frame were Galilean, a Lissajous-Bowditch ellipse. In this case, it will be subjected to a precession of its principal axis, and will have one of the orbits presented in the equations ( 13). Note that it can be interesting to use these initial conditions to improve, in certain cases, the precision of the measurements [10]. The trajectory will then be, in the case of the small angle approximation, one of the hypocrotoid curves, depending on the norm of the coriolis force. We present in Figure 3.3 some of them. Note that the curves are closed if n is a rational number.


Pascal Cachier
email: pascal.cachier@sophia.inria.fr

Xavier Pennec
email: xavier.pennec@sophia.inria.fr

3D Non-Rigid Registration by Gradient Descent on a Gaussian-Windowed Similarity Measure using Convolutions

come

Introduction

In the context of this paper, the goal of non-rigid registration is the retrieval of the motion that has occurred between two volumetric images of the same organ. Because we are working on two images I and J, we have to rely on features of these images. Then, the registration process consists in finding a transformation T that matches well the features, with the hope that these features are significant enough to justify the assumption "features are matched hence homologous physical points are matched." Numerous features are possible, such as points, curves, or intensity. In this paper, we focus on intensity-based methods.

When the intensity is conserved during the mo-tion, an appropriate similarity energy1 is the sum of square differences (SSD). Its gradient is easy to compute and has been used as the guiding force of numerous non-rigid registration techniques [START_REF] Lester | Non-Linear Registration with the Variable Viscosity Fluid Algorithm[END_REF][START_REF] Christensen | Volumetric Transformation of Brain Anatomy[END_REF][START_REF] Bro-Nielsen | Fast Fluid Registration of Medical Images[END_REF][START_REF] Pennec | Understanding the "Demons" Algorithm: 3D Non-Rigid Registration by Gradient Descent[END_REF]. When this assumption is invalid, for example with magnetic resonance images (MRI) corrupted by a nonuniform bias, we have to use more complex criterions, such as the sum of local correlation coefficients (LCC).

In this case, the minimization is usually done by a block matching technique using an exhaustive search [START_REF] Lau | Non-Rigid 3D Image Registration Using Regionally Constrained Matching and the Correlation Ratio[END_REF][START_REF] Maintz | General Multimodal Elastic Registration based on Mutual Information[END_REF][START_REF] Strintzis | Maximum Likelihood Motion Estimation in Ultrasound Image Sequences[END_REF][START_REF] Gee | Bayesian Approach to the Brain Image Matching Problem[END_REF] or a 0th-order minimization strategy [START_REF] Collins | ANIMAL Validation and Applications of nonlinear registration based segmentation[END_REF]. Both the computation of statistics in each hard block and the search of the minimum are computationally intensive for volumetric images. As a consequence, windows are usually sparse and the minimization is done for translations only.

In this paper, we propose to use an isotropic Gaussian window instead of a hard block for the computation of local statistics. This way, we can compute the criterion and its derivatives at each voxel by global convolutions, in a time independent of the extension of the Gaussian window, thus obtaining dense force fields. Section 2 presents the method used to compute local statistics with Gaussian convolutions. Section 3 presents our two-step non-rigid gradient descent based registration algorithm. Finally, we present results of the algorithm on synthetic and real images in Section 4 and we compare different optimization strategies.

Computation of Local Statistics

When the intensity conservation hypothesis that drives numerous mono-modal non-rigid registration algorithms is violated in medical images, this is mainly due to a non-uniform bias that effects the imaging modality. Therefore, the use of local statistics, and particularly the sum of local correlation coefficients, has been increasingly popular. It has even been suggested that it could drive a multi-modal registration [START_REF] Weese | Gray-Value Based Registration of CT and MR Images by Maximization of Local Correlation[END_REF].

Most of the time, these statistics are computed in a window. The image is divided in some number of square windows that may, or may not, overlap. In its simplest form, the block matching scheme consists in moving each window into several positions, computing a local similarity energy inside the window, and keeping the position of the window leading to the best match. However, it is costly to compute the similarity energy for every position of the window. Hence there is a trade-off between precision and computation time. Generally, one chooses fewer windows than voxels. Also, the motion of the window is often restricted to translation, because adding rotation will multiply the number of estimations, especially in 3D images.

Local Gaussian Correlation Coefficients

In the method we propose, we compute local similarity measures around each voxel. The similarity measure we use in this paper is the correlation coefficient.

Let G p be the window function centered around a voxel p obtained by translation of a discrete, normalized, symmetric kernel G, and I and J • T be the two images defined on N n to register with a transformation T . We define the local mean of I, the local correlation between I and J • T , the local variance of I, and the local correlation coefficient between I and J • T (all of them centered around p) respectively by:

Ip = q∈N n Gp[q].I[q] < I, J • T >p = q∈N n Gp[q]. I[q] -Ip J • T [q] -J • T p σ 2 p (I) = < I, I >p CCp (I, J • T ) = < I, J • T >p σp(I) . σp(J • T )
The sum of local correlation coefficients criterion is then defined by:

LCC(I, J • T ) = p∈N n CC p (I, J • T ) (1)
There are multiple choices for the kernel G (hard block, hard sphere, Gaussian, etc.). Here, we have chosen a Gaussian window. The choice of the sum of local correlation coefficients as a registration criterion, and of the discrete, normalized Gaussian as the window function, has been done in the light of the following considerations.

The use of the correlation coefficient enables to compute a distance between two images that is invariant by an affine rescaling of the intensity range on one of the images. It is therefore insensitive to a uniform linear bias.

Because we use local correlation coefficients, the estimation of the bias is done locally. The assumption implicitly done is that the bias should not vary too much inside the window in which the local correlation coefficient is computed. However, the bias is allowed to vary with the position in the image. Hence, we are able to tackle non-uniform biases. We will show later that the LCC criterion (eq. 1) can actually tackle a bias with a variation up to 50% of the intensity range.

The locality of the correlation coefficients is obtained by introducing a Gaussian weight (or Gaussian window) in the formulas. The Gaussian window fulfills two theoretical needs. The first one is the isotropy2 of the window. Hence, if we compute some statistics around a point p, the influence of another point q on the results depends only on its distance to p. The second one is the gradual decrease of the influence of remote points. Remote points should have less influence than points located near the center of the window. This is similar to the implementation of Bajcsy [START_REF] Bajcsy | Multiresolution Elastic Matching[END_REF], where the correlation criterion is not locally weighted but the images are decomposed onto an Hermite orthogonal polynomial basis.

Finally, we compute local correlation coefficients around every voxel, and not only a subset of the voxels. Thus, we do not have any privileged point or location in our images, and our measure is invariant by any integer translation.

Practical Computation using Convolutions

Fortunately, we are able to compute the local statistics and their derivatives at a relatively low computational cost, because the windows are not treated separately but at once, using convolutions. Using the notations previously introduced, local statistics around every voxels can be computed altogether with one or two Gaussian convolutions (denoted by G * ):

Ip = (G * I)[p] σ 2 p (I) = I 2 p -I 2 p = G * (I 2 )[p] -(G * I) 2 [p] < I, J • T >p = I.J • T p -Ip.Jp = G * (I.J • T )[p] -(G * I)[p].(G * J • T )[p]
Note that this can actually be done for any symmetric window function G. However, the Gaussian win-dow, in addition to its theoretical properties presented above, has two practical advantages:

1. The convolution by a Gaussian is separable, i.e., we only have to smooth the image successively with three 1D Gaussian kernels to smooth a 3D image with a 3D Gaussian kernel.

2. The 1D Gaussian convolution can be approximated by a recursive filter, yielding a fast convolution in a time independent of the standard deviation of the Gaussian window [START_REF] Deriche | Recursively Implementing the Gaussian and Its Derivatives[END_REF].

Optimization of the Criterion

To minimize the similarity energy, we compute its derivatives with respect to the displacement field, and use a gradient descent minimization. Different optimization strategies will be presented in the next section. If necessary, we use a multi-resolution scheme to avoid being trapped in local minima.

Because we use Gaussian windows, the derivatives of the criterion can be easily computed using some more Gaussian convolutions. For instance, we put in Table 1 the derivative of the sum of local correlation coefficients criterion with respect to the transformation T [p] of a point p (Eq. 2). We also derived an approximation of the previous formula, where all final convolutions are dropped, in order to be faster (Eq. 3 in Table 1).

In the following, the gradient descent using formula (2) is called Local Correlation Coefficients (LCC). The one using formula (3) is called Simplified Local Correlation Coefficients (SLCC).

We compute a gradient with respect to the displacement of each voxel. Hence, we have a dense vector field. The first advantage of doing so is that we can handle local, subtle deformations, independently of the spread of the Gaussian used to compute local statistics. A second advantage is that the local similarity of the two images is computed using the real image of the window3 and not its rigid or translational approximation (Fig. 1).

Numerical Precision Issues

A mathematical equality may hide numerical problems. For example, consider this equality for the local variance: (I -I p ) 2 p = I 2 p -I 2 p . Whereas the lefthand side will only sum positive and relatively small numbers, the right-hand side will compute two big yet close numbers, I 2 p and I 2 p , and then subtract them: we lose several bits of precision. To quantify the loss Figure 1: On the left, the common block matching scheme assumes that motion is uniform on the entire window and equal to the motion of its center. On the right, using a dense displacement field allows us to use the real image of the window to compute the similarity energy.

of precision, we have compared the results computing the local variances of the image 2(a) using the lefthand side and using the right-hand side. We have used real numbers coded on 4 bytes, and as a preprocessing we have subtracted the mean of the whole image. It appears that the error in the computation of the local variance using the difference of mean relatively to the accurate method has an average of 0.14%±0.83% and a median error of 0.00033%. This means that we have a good numerical precision everywhere, except for strong outliers which appear to lie on the borders. They are due to the conjugate effects of the very low variance on these plain black regions, and the boundary effects of the filter. They are mostly less than 1%, but reach up to 40% in the 4 corner pixels.

In our experiments, these errors never had a meaningful incidence on the results, but we should be aware that it can potentially induce some significant distortions in the deformation field at the image borders.

A Two-Step Registration Algorithm

Non-rigid registration using free form transformations is an ill-posed problem that needs regularization. Within the medical image registration community, regularization is often done by choosing a smoothing energy in addition to the similarity energy, and by minimizing both together. Of course, there is a trade-off between the similarity energy, reflected by the visual quality of the registration, and the smoothing energy, reflected by the regularity of the transformation (the term "regularity" should be taken in its broadest sense, since the smoothing energy may allow occasional discontinuities in the displacement field [START_REF] Hellier | Medical Image Registration with Robust Multigrid Techniques[END_REF]). Therefore, we show the result of a non-rigid registration using the deformed image and the transformation itself.

In the regularization theory framework, one min-

∂LCC(I, J • T ) ∂T [p] = I[p].G * 1 σ q (I).σ q (J • T ) [p] -J • T [p].G * < I, J • T > q σ q (I).σ 3 q (J • T ) [p] + G * < I, J • T > q J • T q σ q (I).σ 3 q (J • T ) - I q σ q (I).σ q (J • T ) [p] ∇J(T [p]) (2)
≈ 1 σ p (I).σ p (J • T ) I[p] -I p -J • T [p] -J • T p < I, J • T > p σ 2 p (J • T ) [p] ∇J(T [p]) (3)
Table 1: Derivative of the sum of local correlation coefficients criterion with respect to the transformation T [p] of a point p (equ. 2) and an approximation without Gaussian convolutions to speed up the computations (equ. 3).

imizes a weighted sum of both energies E 1 + λ.E 2 , E 1 and E 2 being respectively the similarity and the smoothing energies, and λ being a trade-off coefficient. This formulation has proven to be successful for data approximation, and has been used for various approaches of non-rigid registration algorithms [START_REF] Hellier | Medical Image Registration with Robust Multigrid Techniques[END_REF][START_REF] Ferrant | 3D Image Matching using a Finite Element Based Elastic Deformation Model[END_REF].

Another widely used method attempts to separate the image measure from the transformation measure, and could be compared with the approach of game theory. It consists in alternatively decreasing the similarity energy and the smoothing energy. This includes in particular most of block-matching algorithms [START_REF] Lau | Non-Rigid 3D Image Registration Using Regionally Constrained Matching and the Correlation Ratio[END_REF][START_REF] Maintz | General Multimodal Elastic Registration based on Mutual Information[END_REF][START_REF] Strintzis | Maximum Likelihood Motion Estimation in Ultrasound Image Sequences[END_REF][START_REF] Gee | Bayesian Approach to the Brain Image Matching Problem[END_REF][START_REF] Collins | ANIMAL Validation and Applications of nonlinear registration based segmentation[END_REF] and some optical-flow-based techniques [START_REF] Thirion | Image matching as a diffusion process: an analogy with Maxwell's demons[END_REF][START_REF] Horn | Determining Optical Flow[END_REF]. Note that a regularization-theoretic approach may lead in practice to a two-step algorithm [START_REF] Horn | Determining Optical Flow[END_REF].

In this paper, we use the alternated-minimization method. We now detail in turn minimization of the two energies.

Minimization of the Similarity Energy

The image registration step consists in minimizing the similarity energy using a gradient descent [START_REF] Dennis | Numerical Methods for Unconstrained Optimization and Nonlinear Equations[END_REF]. We have tested two strategies: a Gauss-Newton-like scheme, and a linear search.

Gauss-Newton When we want to minimize a nonlinear least squares function z(x) T .z(x) = i (z i (x)) 2 , the Gauss-Newton scheme realizes a second order gradient descent using an approximation of the Hessian matrix from the gradient of this criterion to minimize the computations. The evolution equation is:

x t+1 = x t -(J(x t ) T J(x t ) + λ.Id) -1 J(x t ) T z(x t )
where J = (∂z i /∂x j ) is the Jacobian matrix of z (the matrix of the first order derivatives). In the case of the SSD criterion on two images I and J, this formula adds at each iteration the following vector field to the current displacement field [START_REF] Pennec | Understanding the "Demons" Algorithm: 3D Non-Rigid Registration by Gradient Descent[END_REF]:

u[p] = (I[p] -J • T [p]) ∇J • T [p] 2 + λ ∇J • T [p]
which can also be written:

u[p] = - 2E 1 [p] ∇E 1 [p] 2 + 4λE 1 ∇E 1 [p] (4)
where

E 1 [p] = (I[p] -J • T [p]) 2 is our local similarity energy and ∇E 1 [p] = 2(J • T [p] -I[p])∇J • T [p]
its gradient with respect to the transformation T [p]. To extend this scheme to our non quadratic criterion (the sum of local correlation coefficients), we simply use formula (4) with the similarity energy LCC(I, J • T ) described in Equation 1. Even if the Gauss-Newton scheme was developed for non-linear least squares problems, we will see in the experiments section that our adaptation still performs very well.

Linear search

We have also tried a linear search to compare with our Gauss-Newton adaptation. The search is done along the descent direction given by ( 4), which can be seen as a pre-conditioning of ∇E 1 . The use of a linear search is computationally more expensive since it requires several estimations (typically 8) of the criterion at each iteration.

Minimization of the Smoothing Energy

The regularization of the transformation is done here by decreasing by a certain amount the following stretch energy of the displacement field U = Id -T :

E2(U) = ∂U ∂x 2 = ∂U ∂x 2 + ∂U ∂y 2 + ∂U ∂z 2
This is done by Gaussian smoothing [START_REF] Morel | Variational Methods in Image Segmentation[END_REF], the spread of the kernel growing with the amount of reduction done.

Experiments

To evaluate the different methods previously presented, we first ran a series of experiments on 2D data with synthetic deformations. This allows us to exhibit some key features about each method and criterion while evaluating the accuracy of the non-rigid registration. Then, we show that the algorithms perform similarly on real 3D data. All experiments were performed on a 450 MHz Pentium II and used a Gaussian window with a standard deviation of 4 voxels to compute the local statistics.

Synthetic 2D experiments

In the first experiment, we took a "noise-free" 211×228 MRI slice (obtained by anisotropic diffusion), and we deformed it with a known synthetic transformation meant to be highly non-rigid yet local. Then, we added a white Gaussian noise of standard deviation 3 to both the deformed (Fig. 2(b)) and non-deformed image (Fig. 2(a)).

We tried the Gauss-Newton scheme with the standard SSD, the LCC, and the SLCC formula. We also ran a linear search using the LCC criterion, to compare with the Gauss-Newton-like minimization. 20 iterations were allowed, which is sufficient for the convergence of the Gauss-Newton scheme; the linear search scheme converged around 60 iterations.

In the second experiment, we also added a strong non-uniform bias to the deformed image, ranging linearly from 0 (at the top left corner) to 130 (at the bottom right corner), while maintaining the intensities in the range [0; 255] (Fig. 2(c)).

We have used a Gaussian kernel with a standard deviation of 1.4 for the smoothing. A rough preliminary study has shown that best results were obtained with this standard deviation for these couple of image.

The result of the registration is evaluated by the mean Euclidean distance between the synthetic transformation T * used to deform the image and the transformation T found by the registration algorithm:

D1(T ) = 1 #I p∈I T [p] -T * [p]
where #I is the number of pixels of the image. Since it is impossible to recover the transformation in the black regions outside the head, we have also computed their mean distance inside the head:

D2(T ) = 1 #Ω p∈Ω T [p] -T * [p]|
where Ω is a mask obtained by manual segmentation.

Influence of the minimization scheme First, we observe that the Gauss-Newton-like scheme, which is a special adaptation of a Gauss-Newton for non quadratic criterions, has a proper minimization behavior. The accuracy results are more or less equivalent with the linear search for LCC and SLCC (see Table 2), but the linear search is slower because it needs several computation of the criterion (typically, 8) at each iteration. Table 2: Summary of synthetic experiment results: comparison of accuracy and computation times of algorithms with respect to the criterion and the minimization method used.

Influence of the criterion The sum of local correlation coefficient performs almost as well as the classical SSD on unbiased images. However, as expected, the SSD criterion completely fails for biased images while LCC and SLCC still perform equivalently, showing their robustness to non-uniform bias. The additional computation time due to the use of the local correlation coefficients criterion instead of the SSD criterion is reasonable: a factor 1.75 for SLCC to a factor of 2.50 for LCC. These ratios should actually be smaller when more complex physical models are used for smoothing the deformation field.

Influence of the derivative simplification (SLCC vs LCC) Surprisingly, the algorithm performs slightly better with the derivative simplification (SLCC) than with the exact derivative (LCC), while we were expecting the contrary. We have no explanation for this yet. This result may be due to the fact that we are alternating minimizations and not minimizing a sum of energies. As far as the computation times are concerned, SLCC runs in two thirds of the computation time of LCC.

Analysis of the deformation field errors Errors are mainly due to the inability to totally recover the motion in homogeneous areas and the component of the motion that is tangential to the boundaries of the image and makes them slide. Errors on the component of the motion normal to boundaries are generally much smaller (Fig. 4). Therefore, the result image, whose aspect is somewhat insensitive to the tangential and homogeneous-area errors, looks very much like the target image (Fig. 3).

Real 3D experiment

In this experiment, we have registered two 52×89× 45 MRI centered on the ventricles of a patient who is suffering from focal aphasia. This disease is seen by a dramatic increase of the size of the ventricles between the first MRI 5(a) and the second MRI 5(b) taken a year later. This pair of images presents a bias as well as a contrast change.

The original images were previously rigidly registered with the robust algorithm of [START_REF] Ourselin | Automatic Alignment of Histological Sections[END_REF]. One of the original image does not show the top of the head, which explains why a part of the skull is missing at the bottom of the rigidly registered image 5(a).

We used a Gauss-Newton-like gradient descent with the SLCC formula to register the images, and a Gaussian with a standard deviation of 1.8 to smooth the transformation.

The results presented in Fig. 6 have been obtained in 103 seconds on a 450MHz Pentium II. They show a good match of the structures, as well as a smooth transformation, with deformations localized near the ventricles and the missing part of the skull. This later movement is due to the lack of specific treatment of occlusions in our algorithm. Results using the LCC formula are almost the same, but run in 127 seconds.

A remarkable fact is that an isovalue of the Jacobian of the recovered transformation almost perfectly segments the ventricles which have grown between both acquisitions. The dilatation ratio recovered by the transformation is the same for big structures as well as small structures.

We have also registered these images using the SSD criterion, and a uniform correction of the bias and contrast change, i.e. we have set I = (I -I p ) * σ(J)/σ(I)+J p and registered I and J, I being image 5(a). The registration ran in 64 seconds. The results are poor: this is due to the fact that the uniform intensity correction is not efficient. For example, the background intensity is 0 in J but -9 in I , which explains the problems near the skull. The intensity value of white matter is also underestimated (100 instead of 120). The uniform intensity correction does not work because the histogram is altered by changes (the skull occlusion and mainly the growth of the ventricles) other than the bias and contrast change.

Conclusion

In this paper, we showed how to use the derivatives of the sum of local correlation coefficients (LCC) everywhere in the image to drive a non-rigid registration algorithm. This is computationally feasible thanks to the Gaussian window because we can compute the derivatives very efficiently using separable Gaussian filters. We also compared two gradient de-scent techniques, an adaptation of the Gauss-Newton method and a linear search, to optimize this non quadratic criterion.

We found that the LCC and its approximation SLCC are almost as accurate as the SSD on images without bias. They work almost as well with biased images, where the SSD fails, showing the robustness of this criterion with respect to non-uniform bias. The LCC and the SLCC methods demonstrate equivalent accuracies, but SLCC is faster than LCC: they were respectively 1.75 and 2.50 times slower than the SSD in our experiments.

A few points still need to be investigated to fully control and automate the algorithm. Firstly, we have to study the impact of the boundary conditions of the filters on the computations of the similarity measure and the gradient. The spread of the Gaussian window used to compute the local statistics could be chosen automatically. Last but not least, we should study and compare the behavior of our technique using a minimization of the sum of the similarity and the smoothing energies.

We believe that numerous non-rigid registration algorithms that use dense (i.e., at each voxel) forces to guide the registration could be extended from the SSD to the LCC criterion using the approach presented here. Moreover, even if we presented the formulas and the experiments with the sum of local correlation coefficients, the method is intrinsically generic and could be applied to many other local criterions, the only limitation being the ability to compute the derivative of the criterion.

Figure 2 :

 2 Figure 2: Synthetic images and transformation. From left to right: Image 2(a) has been transformed into image 2(b) using a deformation depicted on image 2(d), and white Gaussian noise (σ = 3) has been added afterwards. Image 2(c) has been obtained by adding a bias varying linearly from 0 at the top left corner to 130 at the bottom right corner to image 2(c) (the original image intensity ranges from 0 to 250).

Figure 3 :

 3 Figure 3: Result of non-rigid registration of image 2(a) toward image 2(c) using the SLCC and Gauss-Newton strategy. Left image should be compared to the target image 2(c). Structures are almost perfectly matched. Right image should be compared to the synthetic transformation 2(d). The transformation is recovered in regions carrying information.

Figure 4 :

 4 Figure 4: Errors on the displacement field for the region delimited on the left image. Middle image shows the errors along the x-component, normal to the boundary of the skull. Errors are around 0.1 pixel. Right image shows the errors along the y-component, tangential to the boundary of the skull. Local maxima of the error are between 2.6 and 3.0 pixel (local maxima of the initial displacement is 3.0 pixel). The tangential component of the displacement is not recovered and constitutes the main source of errors.

Figure 5 :

 5 Figure 5: Rigidly-registered 3D MR images of the ventricles. Notice the bias and the deformation in relation with focal aphasia.

Figure 6 :

 6 Figure 6: Result of non-rigid registration of the images 5 using the SLCC+Gauss-Newton algorithm. Left Image shows the image result, which should be compared to the image 5(b). Ventricles are good, but grey matter is difficult to match, and the skull occlusion perturbs the deformation. Right Image shows the transformation result. It is smooth, and the biggest deformations are inside the ventricles, and near the missing part of the skull.

Figure 7 :

 7 Figure 7: An isovalue of the recovered transformation segments the ventricles.

Figure 8 :

 8 Figure 8: Results using the SSD + Gauss-Newton algorithm with uniform intensity correction. Due to the growth of the ventricles, the uniform intensity correction underestimates the intensities. This is enough for the SSD to fail.

For convenience, we shall talk about similarity energies instead of similarity measures. This terminology supposes that the energy should be positive and takes zero as a minimum. A similarity measure S with an upper bound S * can be transformed into a similarity energy E by setting E = S * -S.

The Gaussian should be isotropic in the physical space. If the voxels themselves are not isotropic, then the Gaussian window follows this anisotropy in the image space.

If we use a recursive implementation for the Gaussian, the window actually fills the entire image.

Acknowledgments: We thank S. Ourselin for his proofreading and fruitful discussions.


Monireh Dehkhodaeia

Mehdi Sahihia

Hadi Amiri Rudbaria

Dr Mehdi Sahihi
email: m.sahihi@chem.ui.ac.ir

Spectroscopic and molecular docking studies on the interaction of Pd(II) & Co(II) Schiff base complexes with β-lactoglobulin as a carrier protein

Keywords: Schiff base complex, β-lactoglobulin, interaction of BLG, Molecular docking

Introduction

The milk is one of the most popular and of course the most essential beverage. This food source is rich of proteins, vitamins and the other necessary ions for body healthy. The protein part of milk is contained of casein (78.3%), whey proteins (19%) and the other proteins (2.7%). Whey proteins are significantly, used in the food and pharmaceutical industry due to their high nutritional value and excellent functional properties. The most important and abundant part of whey proteins is beta-lactoglobulin (BLG) with molecular weight of 18.400 Da [START_REF] Papiz | The structure of b-lactoglobulin and its similarity to plasma retinol-binding protein[END_REF]. This globular protein is contained of 162 amino acids that they are folded in structure including one α-helix, nine antiparallel β-strand with eight β-sheetforming a cone barrel shape with a hydrophobic pocket. This hydrophobic pocket is known as calyx and it is a very convenient location for binding ligands and drugs with very low solubility in water. Generally, BLG has three binding sites for suchlike ligands. Results of previous studies show that vast range of hydrophobic molecules such fatty acids, cholesterol, porphyrin, bilirubin, aromatic hydrocarbons, vitamin D3, anticancer and chemotherapy drugs can bind to BLG for transporting in body.Also, BLG is known as a carrier protein with controlled release ability [START_REF] Reddy | Structural and conformational basis of the resistance of. beta.-lactoglobulin to peptic and chymotryptic digestion[END_REF]. On the other hand, finding novel compounds with therapeutic properties, especially anticancer compounds, for binding to carrier proteins has biological and pharmacological importance [START_REF] Moghaddam | Comparison of the binding behavior of FCCP with HSA and HTF as determined by spectroscopic and molecular modeling techniques[END_REF][START_REF] Omidvar | Investigations with spectroscopy, zeta potential and molecular modeling of the non-cooperative behaviour between cyclophosphamide hydrochloride and aspirin upon interaction with human serum albumin: Binary and ternary systems from the view point of multi-drug therapy[END_REF][START_REF] Omidvar | Studies on the antagonistic behavior between cyclophosphamide hydrochloride and aspirin with human serum albumin: time-resolved fluorescence spectroscopy and isothermal titration calorimetry[END_REF][START_REF] Hosainzadeh | Probing the interaction of human serum albumin with bilirubin in the presence of aspirin by multi-spectroscopic, molecular modeling and zeta potential techniques: Insight on binary and ternary systems[END_REF].Recently, attention to synthesis of metal based drugs has been enhanced. Schiff-bases are compounds with imine (C=N) functional group, which are formed from condensation reaction between aldehyde and amine. Researchers have focused on these compounds due to their biological aspects like antifungal, antitumor, antiviral, antioxidant and so on [START_REF] Da Silva | Schiff bases: A short review of their antimicrobial activities[END_REF]. Biochemists believethat the existence of imine group in Schiff-basesligands is responsible for these biological effects. Also, the metal ions play important role in pharmaceutical properties of these low water soluble compounds, e.g.importance of Co(II) in vitamin B12 or Pd(II) role in anticancer and antitumor agents. With consideration the above-mentioned points, studying on the interaction of Schiff base metal complexes and carrier proteins is very attractive for researches. In this study we investigated the interaction of PdL2and CoL2Schiff-base complexes which are synthesized from condensation reaction of salicylaldehyde and tert-butylamine (Scheme S1), (these complexes are prepared and characterized in our research group recently [START_REF] Khorshidifard | Cobalt (II), copper (II), zinc (II) and palladium (II) Schiff base complexes: Synthesis, characterization and catalytic performance in selective oxidation of sulfides using hydrogen peroxide under solvent-free conditions[END_REF]) with BLG using UV-Vis absorption, fluorescence spectroscopic and molecular docking methods.

Experimental.

Material

All of the materials for synthesis of the complexes were purchased from Merck-Co and were used without further purification. Bovine BLG type A was purchased from Sigma-Aldrich. All salts used for buffer preparation were analytical grade and were dissolved in double distilled water. All of the solutions were used freshly after preparation.

Preparation of compounds and protein stock solutions.

For optical measurements, all solutions were prepared in phosphate buffer (50 mM) at pH=7. A stock solution of BLG was prepared by dissolving the desired amount of BLG in phosphate buffer.The concentration of BLG was determined using the value of 17,600 M-1.cm-1 for the molar absorption coefficient (ɛ) at 279 nm.The solutions of the compounds were first prepared in methanol and DMSO for CoL2 and PdL2, as co-solvents, respectively. And then diluted suitably with phosphate buffer to required concentrations. The volume of co-solventsnever exceeded 3% (v/v), so their effectswere negligible. All the solutions were used freshly after preparation.

Fluorescence spectroscopy study

The fluorescence spectra were obtained at room temperature on a Shimadzu RF-5000 spectrofluorometer in quartz cell. For studying thefluorescence quenching of BLG by adding PdL2 and CoL2complex solutions, 2mL of BLG solution (10-5 M) was placed in the cell and various amounts of the compounds (2×10-4 M) were added to the cell.In each measurement, the mixture was allowed to incubate for 2 min. The fluorescence intensity was measured with excitation wavelength at 295 nm and emission wavelength range of 300-450 nm. Moreover, the measured fluorescence intensities were corrected for the dilution and the inner-filter effects. To eliminate the inner filter effect, which is happen in high concentrations of protein, absorption measurements were carried out at the fluorescence excitation and emission wavelengths. The extent of this effect can be roughly evaluated with the following relationship:

Where, Fcorr and Fobs are the corrected and observed fluorescence intensities, respectively, and Aex and Aem are the absorption of compoundat excitation and emission wavelengths, respectively.

UV-Vis absorption spectroscopy study

The UV-Vis absorption spectroscopy measurements are carried out to more investigate BLG-binding of the compounds at room temperature. Electronic absorption spectraare recorded on a PERKIN ELLMER Lambda 265 spectrophotometer. For obtaining the UV-Vis absorption spectra 2mL of BLG solution (2×10-5 M) was placed in the cell and various amounts of the compounds (2×10-4 M) were added and spectrum changeswas followed in the absence and presence of the compounds on 279 nm. All intensities were corrected for dilution effect in fluorescence andUV absorption experiments.

Molecular docking procedure

Molecular docking is valuable method to study the interaction of protein and ligand and also, is one of the original tools in computational design of drugs. In this work, we docked the complexes with BLG as an important carrier protein. The 3D structures of the metal complexes were obtained using the CIF files of their X-ray crystal structures [START_REF] Khorshidifard | Cobalt (II), copper (II), zinc (II) and palladium (II) Schiff base complexes: Synthesis, characterization and catalytic performance in selective oxidation of sulfides using hydrogen peroxide under solvent-free conditions[END_REF]. The CIF files were converted to the PDB format by using the Mercury software (http://www.ccdc.cam.ac.uk/). The crystal structure of BLG (PDB ID:1CJ5) was taken from the Brookhaven Protein Data Bank (http://www.rcsb.org/pdb). The resolution of this file was2.28 Å. Water molecules of the .pdb fileswere deleted, missing hydrogen atoms and Gasteiger charges were added. Flexible-ligand docking was carried out by AutoDock 4.2.5.1 molecular-docking program using the implemented empirical free energy function and the Lamarckian Genetic Algorithm. The Gasteiger charges were added to prepare the macromolecule input file for docking and the Auto Grid was used to calculate Grid. For the docking of metal complexes with BLG, a blind docking with 126 lattice points along X, Y, and Z axes was performed to find the binding site of complexes on BLG with a grid point spacing of 0.375 Å, to allow the complex to rotate freely. In the next step, the centre of the grid box was located at the binding site and the second docking was performed using a cubic box with 60×60×60 dimensions. 250docking runs with 25,000,000 energy evaluations for each run were performed.

Result & discussion

Fluorescence quenching analysis

Whereas, BLG known as a carrier protein, understanding of its physicochemical and pharmaceutical features are important steps in investigation of its biological activity. This major bovine whey protein contains remarkable properties that have been previously pointed out. Moreover, its thermodynamic stability in the presence of some surfactants has been investigated, recently. [START_REF] Marouzi | Study on effect of lomefloxacin on human holo-transferrin in the presence of essential and nonessential amino acids: Spectroscopic and molecular modeling approaches[END_REF][START_REF] Ahmadi | A comparison study of the interaction between βlactoglobulin and retinol at two different conditions: spectroscopic and molecular modeling approaches[END_REF]Chamani, et al., 2006;[START_REF] Chamani | Comparison of the conformational stability of the non-native αhelical intermediate of thiol-modified β-lactoglobulin upon interaction with sodium n-alkyl sulfates at two different pH[END_REF][START_REF] Chavoshpour-Natanzi | Structural stability of βlactoglobulin in the presence of cetylpyridinum bromide: spectroscopic and molecular docking studies[END_REF]. Moreover, the appearance of isosbestic points signifies that equilibrium is established between BLG and the compounds and reflects the formation of the ground state complex between them. Whereas, estimating the binding constant of ligand to protein is an important parameter in investigation of interaction of, itcould be given by the following equation that it is known as Scatchard equation:

Where F0, F and [Q] are the same as the parameters in the Stern-Volmer equation. Kb is binding constant that is obtained from exponential y-intercept. In this equation (n) is the number of binding site of the compound on the protein that it is calculated from the plot of ln((𝐹0 -𝐹)/𝐹) versus ln[Q] as plot slop. This plots at T=298 K are presented in Figures1b and 1e. Also, these plots for other temperatures are presented in Figures S2 andS3. The estimated values of Kb, n, Ksv and kq are listed in Table1. The Stern-Volmer relationship, allows us to explore the kinetics of a photophysical intermolecular deactivation process. However, the binding constant, or association constant, is a special case of the equilibrium constantK. It is associated with the binding and unbinding reaction of receptor and ligand molecules and allows us to explore the thermodynamic of binding process. Hence, the quantities of these constants are different. Also, fluorescence spectroscopy is more sensitive than the absorption spectroscopy to assess the binding of drugs to proteins. Hence, the obtained binding constants using fluorescence titration method are more accurate. Generally, binding constant should not be too high which drug cannot separate hereinafter binding to carrier protein and it should not be too low which cannot bind well to carrier protein and achieve to destination. BLG-binding constants for both of the compounds are in a good range (2-6×103). Also, in this study, the values of kq were obtained about 1010 M-1.S-1at room temperature.This observation supports that fluorescence quenching of BLG occurs by static mechanism. The charts of 𝐹0/𝐹 vs [Comp] for both of the compounds in 298 K are present in Figure S1.

Analysis of UV-Vis absorption of BLG

UV-Vis absorption spectroscopy is another technique to investigate the interaction of ligand-protein systems. Although sensitivity of this method is not comparable with fluorescence spectroscopy, but it can be employed as a supplement technique. Concentrations of compounds and BLGarechosen somehow, to cover mole ratio the ([Comp]/[BLG])= 0-5. By using Benesi-Hildebrand plot and drawing the plot of 1/(A-A0) versus 1/[M],binding constant can be estimated as ratio of y-intercept to slop:

Where A0 and A are the absorbance of BLG in the absence and presence of the compounds, respectively. Amax is the obtained absorbance at saturation and [M] is the concentration of the compound. The estimated values of Kb for PdL2and CoL2 complexes at T=298 K are 1.73×104 and 1.72×104, respectively. The UV-Vis spectra of BLG in the presence of different amounts of PdL2 and CoL2 complexes are shown in Figures 2a and2c.

Energy transfer between BLGand complexes.

Fluorescence resonance energy transfer (FRET) or Förster's theory is a non-destructive spectroscopic method that can estimate distance from the bound ligand to the fluorophore groups of biomacromolecule. The energy transfer could take place through direct electrodynamic interaction between the primarily excited molecule and its neighbours. The following conditions determine the occurrence of energy transfer between a donor and an acceptor: (i) the relative orientation of the donor and acceptor electric dipole moments; (ii) the overlap between the emission spectrum of the donor and the UV-Vis spectrum of the acceptor; and (iii) the distance between the donor and the acceptor [START_REF] Bai | Fluorescence resonance energy transfer between bovine serum albumin and fluoresceinamine[END_REF]. According to the Förster's dipole-dipole non-radiative energy transfer theory, the efficiency of energy transfer (E) is related not only to the distance between the acceptor and the donor, but also to the critical energy transfer distance (R0), as described by the following equation:

Where F and F0 are the fluorescence intensity of the donor in the absence and presence of equal amount of acceptor, r is the distance from the bound ligand to the tryptophan residue, and R0 is the Förster critical distance at which 50% of the excitation energy is transferred to the acceptor and can be calculated by the following equation:

Where,K2 is factor describing the relative orientation of the transition dipoles of the donor and acceptor (for a random orientation as in fluid, K2 = 2/3), N is the average refractive index of the medium in the wavelength range where spectral overlap is significant (N is 1.3598), 𝜑 isthe fluorescence quantum yield of the donor (φ is 0.08), and overlap integral 𝐽 expresses the extent of overlap between the normalized fluorescence emission spectrum of the donor and the acceptor absorption spectrum. 𝐽is given by the following equation:

Where, F(λ) is the fluorescence intensity of the fluorescent donor at wavelength λ and is dimensionless, ε(λ) is the molar absorption coefficient of the acceptor at wavelength λ. As it is shown in Figures 1c and1f, the overlap between the fluorescence emission of BLG and the absorption spectra of the compounds is considerable. The calculated values for E, r and R0 parameters are collected in Table 1.

Analysis of thermodynamic parameters

Results of the reports in literature show that there are four non-covalent effective forces between compounds and biomacromolecules, i.e. hydrophobic, H-bond, Vander Walls and electrostatic interactions. According to the Van't Hoff equation, thermodynamic parameters are estimated by the following equation:

In both of these equations, R is gas public constant, T is the experimental temperature, K is binding constant at the corresponding temperature, ΔG°, ΔH° and ΔS° are the standard Gibbs free energy, the enthalpy and the entropy changes, respectively. ΔH° and ΔS° are estimated from the plot lnK vs 1/T as slop and y-intercept, respectively. These plots for CoL2 and PdL2 complexes are present in Figures S4a andS4b, respectively. Also ΔG° is obtained by using equation ( 7). According to the Rass and Subramanian report [START_REF] Ross | Thermodynamics of protein association reactions: forces contributing to stability[END_REF]:

The obtained results in Table 2show that, both of ΔH° and ΔS°parameters are positive.Hence, hydrophobic interactions have main role in the stability of BLGcomplexes systems.Also, negative values of ΔG°at alltemperatures indicates that the interaction of BLG and PdL2 or CoL2 is spontaneous.

Molecular docking study

Studies on the interaction between BLG and various compounds show that there are three potential binding sites on BLG: the internal cavity of the β-barrel that is known as main binding site, the surface hydrophobic pocket in a groove between α-helix and βbarrel and the outer surface near TRP19-ARG124. Figures 3a and3b indicate the binding site position for PdL2and CoL2with their chosen clusters, respectively. The obtained results showed that PdL2 complex has located in secondary binding site of BLG, between α-helix and β-barrelwith standard binding free energy (ΔG°) of -6.3 kcal.mol-1 (Kb= 4.17×104 M-1). This score is related to the docking cluster with minimum ΔG° (see Figure 3a).Also, three π-cation interactions have been formed between N and H atoms of LYS138 residue with C atoms of phenyl ring of PdL2 complex. The residues around of this complexare GLN5, THR4, VAL3, TYR2, LEU95, LEU104, LYS138, TYR102, ALA139, ALA142, LEU140 and LEU143. Also, intermolecular energy includingH-bonding interactions, van der walls, desolvation and electrostatic energies,was estimated as -6.85 kcal.mol-1.Figure 3b showsthe position of CoL2 complexon BLG.It isplaced in entrance of the calyx and is stabilized by one hydrogen bond with SER116 residue of BLG. The amino acids with major role in the stability of CoL2-BLG system are GLU112, GLN115, PRO113, ASN109, SER116, ASN88, LEU31, PRO38 and LEU39. The intermolecular energy and standard binding free energy (ΔG°) for CoL2 complex are -6.42 and -5.87 kcal.mol-1 (Kb= 2.02×104 M-1), respectively. These energies are related to the docking cluster with minimum energy (see Figure 3b). Supportingly, differences between the obtained binding constants using fluorescence spectroscopy and molecular docking methods are related to their different natures. However, PdL2 shows more affinity than CoL2 for binding to BLG in all of the experimental and theoretical methods. Furthermore, the distances between PdL2 and CoL2 with Trp 19 residue are 1.21 and 1.63 nm, respectively. These obtained data are in consistent with obtained efficiency of energy transfer (E) using FRET experiments (Table 1).

Conclusion

Two metalSchiff-base complexes(CoL2 and PdL2)have been investigated to obtain binding affinity to BLG using UV-Vis, fluorescence spectroscopy and molecular docking methods. The obtained results of the experimental section showed that the compounds bind to BLG as PdL2˃CoL2. Also, to calculate thermodynamic parameters, fluorescence spectroscopy was done at three temperatures (298, 308 and 318 K).The obtained results indicated negative values for ΔG° and positive values for ΔS° and ΔH°,at all temperatures. Also, molecular docking studies revealed that hydrogen bond,hydrophobic and π-cation interactions have the important role in the interaction of these compoundswith BLG.

  Fluorescence spectroscopy is one of the trusty methods to check the interaction of ligand and protein. This high sensitive technique can indicate accessibility of quencher to fluorophoresof protein. Fluorescence quenching could proceed through different mechanisms, usually classified as dynamic and static quenching. Static quenching mode refers to the formation of a non-fluorescent ground state complex, whereas, dynamic quenching mode is resulted from the collision between the fluorophore and quencher in excited state. Both static and dynamic processes are described by the Stern-Volmer equation as follow: In this equation, F0 and F are the fluorescence intensities of BLG in the absence and presence of the compounds. Ksv is the Stern-Volmer quenching constant, kq is the quenching rate constant of BLG, τ is the average lifetime of BLG without quencher which is typically equal to 10-8 s for biomacromolecules and [Q] is the concentration of quencher that quencher is the synthesized compounds here. Ksv is determined from the plot of 𝐹0/𝐹 vs [Q]. The fluorescence spectra of BLG in the presence of different amounts of PdL2 and CoL2 complexes at T=298 K were shown in Figures1a and 1d.Also these spectra at T=308 and 318 K are shown inFigures S2 and S3. The fluorescence intensity of BLG gradually decreased by adding different volumes (μL) of compounds in each step.

Figure 1 .

 1 Figure 1. a and d)The quenching effect of schiff-base complexes on BLG fluorescence intensity at 25° C for CoL2 and PdL2, respectively. [BLG]=10-5 M, [Comp]=2×10-4 Mand λex=295 nm.b and e) Plots of ln((F0-F)/F) versus ln[Comp] for CoL2 and PdL2, respectively.c and f) Spectral overlap of the BLG fluorescence with the absorption spectraof CoL2 and PdL2, respectively.

Figure 2 .

 2 Figure 2.a and c) Changing of UV-Vis spectrum of BLG upon its titration with various amounts of CoL2 and PdL2, respectively. [BLG]=2×10-5 M, [Comp]=2×10-4 M.b and d) Plots of 1/A-A0 versus (1/[Comp]) for CoL2 and PdL2, respectively.

Figure 3

 3 Figure 3.a) View of the BLG binding site for PdL2andits selected cluster. b) View of the BLG binding site for CoL2 andits selected cluster. (Small spheres and yellow cones show H-bond interaction and π-cation interactions, respectively.)

Figure

  Figure S2.a and c) The quenching effect of schiff-base complexes on BLG fluorescence intensity fora) CoL2, c)PdL2, respectively,at 35° C. b and d) Plots of ln((F0-F)/F) versus ln[Comp] for CoL2 and PdL2, respectively,at 35° C.

Figure

  Figure S3.a and c) The quenching effect of schiff-base complexes on BLG fluorescence intensity fora) CoL2, c)PdL2, respectively,at 45° C.b and d) Plots of ln((F0-F)/F) versus ln[Comp] for CoL2 and PdL2, respectively,at 45° C.