Cosmological solutions in the Brans-Dicke theory
via invariants of symmetry groups

E. Ahmadi<EMAIL_ADDRESS>K.
<EMAIL_ADDRESS>A<EMAIL_ADDRESS>
Department of Physics, Faculty of Basic Sciences,
Azarbaijan Shahid Madani University, 53714-161, Tabriz, Iran

We proceed to obtain an exact analytical solution of the Brans-Dicke (BD)
equations for the spatially flat ($k=0$) Friedmann-Lamaitre-Robertson-Walker
(FLRW) cosmological model in both cases of the absence and presence of the
cosmological constant. The solution method that we use to solve the field
equations of the BD equations is called the “invariants of symmetry groups
method” (ISG-method). This method is based on the extended Prelle-Singer (PS)
method and it employs the Lie point symmetry, $\lambda$-symmetry, and Darboux
polynomials (DPs). Indeed, the ISG-method tries to provide two independent
first-order invariants associated to the one-parameter Lie groups of
transformations keeping ordinary differential equations (ODEs) invariant, as
solutions. It should be noted for integrable ODEs, the ISG-method guarantees
the extraction of these two invariants. In this work for the BD equations in
FLRW cosmological model, we find the Lie point symmetries,
$\lambda$-symmetries and DPs, and obtain the basic quantities of the extended
PS method (which are the null forms and the integrating factors). By making
use of the extended PS method we find two independent first-order invariants,
in such a way appropriate cosmological solutions from solving these invariants
as a system of algebraic equations are simultaneously obtained. These
solutions are wealthy so that they include many known special solutions, such
as O’Hanlon-Tupper vacuum solutions, Nariai’s solutions, Brans-Dicke dust
solutions, inflationary solutions, and etc.
---

Keywords: Cosmological Solutions, Brans-Dicke Theory, Invariants of Symmetry
Groups Method

###### Contents

  1. 1  Introduction
  2. 2 A short review of the ISG-method
  3. 3 FLRW cosmological model in BD theory
    1. 3.1 A short review of the theory of BD
    2. 3.2 FLRW cosmological model in BD theory
  4. 4 Solving the FBD-equation
  5. 5  Completing the solutions
  6. 6 Special solutions of the FBD-equation
  7. 7 Conclusions

##  1  Introduction

Mach’s ideas strongly inspired Albert Einstein, so he attempted to include
this principle in his theory when formulated the theory of general relativity
(GR). During the years 1907-1916 when Einstein was developing his theory of
GR, he strongly influenced by Mach’s philosophy and tried to include this
principle in his theory in every possible. Finally, he could reflect the
relationship between the geometry of space-time and also the distribution of
matter in the GR equations. Later, it was shown that the Einstein’s field
equations (EFEs) can be include solutions in the absence of matter and in
empty space (Minkowski space-time, Taub-NUT space-time, etc.) [1]. Therefore,
contrary to Mach’s opinion (the stronger form of Mach’s principle), the GR
theory includes the weaker form of Mach’s principle, and hence space-time has
an existence independent of matter. Among physicists, there is a group that
accepts Mach’s opinion and believes that GR, which does not include the
stronger form of Mach’s principle, should be modified to include it. Since the
70s, efforts have been made in this direction, among which we can mention the
scalar-tensor theories (STTs) [2]. The Brans-Dicke (BD) theory of gravity
(sometimes called the Jordan-Brans-Dicke scalar-tenser theory of gravity) is
the simplest version of the STT, which was developed by Pascual Jordan in the
1950s and later by Brans and Dicke in the early 1960s [3, 4]. Actually,
theories about this had already been proposed years earlier by Jordan, Fierz
and Thiery. One of the motivations of these authors was to include the
stronger form of the Mach’s principle in gravity. This principle states that
the inertial mass of a body is due to all the rest of the matter in the
universe. For the first time, the idea of varying the Newton’s gravitational
constant $G$ with the age of the universe was proposed by Dirac [5]. His
purpose in doing this was to explain the weakness of the gravitational force
among the fundamental forces in nature. In the Jordan gravity theory, the
variation of the constant $G$ with space-time was described by a simple scalar
field. However, scientists such as Fierz and Bondi criticized this theory,
because energy was not conserved in it. In 1961, Brans and Dicke proposed [3]
a more complete theory of gravity incorporating a varying gravitational
constant (see also [4]). In the BD theory, the constant $G$ has been replaced
by the inverse of a scalar field $\varphi$. The scalar field $\varphi$ plays
the role of the $G$, and Mach’s principle is revived through it. Also, in the
BD theory a dimensionless parameter $\omega$, which measures the coupling of
the $\varphi$ to matter, is included. In the GR, there is only one tensor
field $g_{\mu\nu}$ describing the geometry of space-time, while in the BD
theory, in addition to the tensor field $g_{\mu\nu}$, a scalar field $\varphi$
is incorporated for reviving Mach’s principle. In 1980, Batista [6] extended
the cosmological model obtained by Brans and Dicke to non-homogeneous models.
He considered density $\rho$ and pressure $p$ as functions of the variables
$r,t$ and showed that for the dust-dominant universe, the same solutions for
$a(t)$ and $\varphi$ are obtained as previously found by Brans and Dicke. In
other words, the assumption of inhomogeneity of dust distribution goes not
change the BD solutions. McIntosh [7] showed that with the spatially flat FLRW
metric, the scalar field $\varphi$, the pressure $p$ and its density $\rho$,
there are solutions for the BD equations. He showed that for solutions of the
BD field equations, it must not be assumed that the $\varphi$, $p$ and $\rho$
are only functions of the variables on which the metric $g_{\mu\nu}$ depend.
This true not just in the cosmological case but in vacuum and other cases as
well. The vacuum solutions of BD equations for the spatially flat FLRW metric
were obtained by O’Hanlon and Tupper [8], Chauvet [9], Cervero [10] and
Lorentz-Petzold [11]. Chauvet also obtained the solutions of the BD equations
parametrically in the case $k=\pm 1$. Uehara and Kim [12] studied the BD
theory in the presence of the cosmological constant $\Lambda$, in such a way
that they obtained solutions for the BD equations in the presence of the
cosmological constant for the spatially flat FLRW metric. They also obtained
an exponential solution for the radiation universe $(p=\rho/3)$ when
$\omega=-3/2$.

Cosmological dynamical equations (DEs) of FLRW model in the framework of GR
theory and many other models in the framework of extended gravity theories
such as BD theory [13], $f(T)$ gravity [14], $f(R)$ theory [15], Rastall
theory [16], etc., have been solved analytically and numerically by different
methods. But, there is no unique solution method that is general, so that it
can be used not only to solve the field equations of BD theory, but also the
field equations of any extended gravitational theories. Note that the
cosmological DEs are coupled non-linear differential equations, which require
special methods. The PS-method, which was later modified and completed by
Chandrasekar et al. [17], is one of these methods. This solution method has
made great progress in the field of theory of non-linear differential
equations. The ability that this method has shown in solving a number of well-
known non-linear differential equations, such as the $2$-dimensional Kepler
problem, the modified Emden equation, and the Helmholtz oscillator [18], is
indicative of the fact that the modified extended PS-method, when combined
with other symmetry methods, it can be a powerful algorithmic solution for
solving non-linear DEs in physics, such as cosmology [19, 20]. We may use this
method to extract independent invariants from the cosmological DEs in a mini-
super space, and hence we have actually solved the problem, because these
invariants can be considered as a system of algebraic equations. Solving this
system of algebraic equations is much easier than solving the system of
cosmological DEs (which are in the form of non-linear differential equations).
It should be noted that the extended PS-method and combined with other
symmetry methods that we call the ISG-method, they are not included in other
branches of science, especially gravity and cosmology. The remarkable point in
this method is that it is difficult to calculate its basic quantities, i.e.,
null form $S$ and integrating factor $R$ in the extended PS-method are related
to the fundamental quantities of other symmetry methods such as the
infinitesimals $\tau$ and $\xi$ in the Lie point symmetry and the function
$\lambda$ in the $\lambda$-symmetries, and DPs in the DPs method [19, 20, 21,
22, 23, 24, 25, 26]. Recently, we have presented the ISG-method in order to
solve the field equations of the spatially flat FLRW cosmological model in the
presence of $\Lambda$ in the context of GR [27]. Here, we employ this method
to obtain the solutions of the field equations of BD theory for spatially flat
FLRW cosmological model in both cases of the absence and presence of the
cosmological constant. We show that when the coupling constant $\omega$ tends
to infinity, the cosmological solutions of GR theory can be obtained from the
analytical solutions of BDEs. For the dust dominated universe, when the
cosmological constant is present, we conclude that the solutions of Uehara and
Kim [20] are special cases of the solutions of BDEs. In addition, we show that
the case $\Lambda=0$ of the solutions of BDEs include Nariai’s solutions, BD
dust solutions, radiation solutions, inflationary solutions, O’Hanlon-Tupper
vacuum solutions and finally general relativity solutions.

The paper is organized as follows. In Sec. 2, we briefly describe the ISG-
method applicable for second-order ODEs in physics, especially in cosmology.
We give the main points of the ISG-method in this section. In Sec. 3, we first
introduce the BD theory of gravity, and then extract the DEs of the FLRW
cosmological model in this theory. Then, by combination these equations, we
will reduce them to a second-order ODE called the Friedmann-Brans-Dicke
equation (FBD-equation). In Sec. 4, we will apply the ISG-method to solve the
FBD-equation. This solution enables us to obtain the cosmological solutions of
the BD theory. In Sec. 5, the solutions of the BDEs are completed so that they
are suitable for cosmological applications. In Sec. 6, we will examine the
cosmological solutions in special cases. Finally, conclusions are reported in
Sec. 7.

##  2 A short review of the ISG-method

In order to introduce the ISG-method [27], let us consider a particle of unit
mass moveing in one-dimensional configuration space $\mathbb{Q}=(q)$. A force
${\bf F}=\phi(t,q,\dot{q}){\partial}/{\partial q}$ is applied to the particle
and gives it an acceleration ${\bf a}=\ddot{q}{\partial}/{\partial q}$.
According to Newton’s second law of motion, the DE of the particle is as
follows:

$\ddot{q}=\phi(t,q,\dot{q}),$ (2.1)

where $\phi$ as a function of $t$, $q$ and $\dot{q}$ is called a force
function of the particle. Here, “.” denotes the total derivative operator with
respect to the time $t$, which is defined as

${}^{.}=\frac{d}{dt}=\frac{\partial}{\partial
t}+\dot{q}\frac{\partial}{\partial
q}+\phi(t,q,\dot{q})\frac{\partial}{\partial\dot{q}}.$ (2.2)

In the force function $\phi(t,q,\dot{q})$, time $t$ is the independent
variable, and coordinate $q$ is the dependent variable. Many physical
phenomena have a DE in the form of equation (2.1). That is, the dynamics of
such phenomena can be imagined as the dynamic of a particle with a unit mass,
which is subject to force ${\bf F}$. The ISG-method provides us an algorithmic
and systematic solution processes for solving second-order ODEs of the form
equation (2.1), in which the force function $\phi(t,q,\dot{q})=P/Q$ is a
fractional function of polynomials $P$ and $Q$ of the variables $t$, $q$, and
$\dot{q}$ with coefficients in the set of complex numbers. In this work, we
use the ISG-method to solve second-order non-linear ODEs like (2.1), in which
the force function is as in those of (2.1). This method can be developed to
solve non-linear ODEs of higher orders and even for system of ODEs of any
order, which are not our discussion here. The general strategy of the ISG-
method for solving second-order ODE (2.1) is to systematically and
algorithmically extract a certain number (equal to the order of the DE, which
is two here) of the independent “first integral” (or first-order “invariant”)
such as $I_{1}(t,q,\dot{q})=c_{1}$ and $I_{2}(t,q,\dot{q})=c_{2}$ associated
to the DE (2.1), where $I_{1}$ and $I_{2}$ are known functions of the
variables $t$, $q$ and $\dot{q}$ and $c_{i}$’s ($i$=1,2) are constants on the
solutions of the DE (2.1). When we obtain these two invariants, the solution
of the DE (2.1) is actually prepared, because one can look at the first-order
ODEs, $I_{1}(t,q,\dot{q})=c_{1}$ and $I_{2}(t,q,\dot{q})=c_{2}$ as a system of
algebraic equations in terms of variables $t,q$ and constants of motion
$c_{1}$ and $c_{2}$. It is easy to solve this system of equations. To solve
this system of algebraic equations, it is enough to remove the variable
$\dot{q}$ between them. For this purpose, one can obtain from one of them, for
example, $I_{1}(t,q,\dot{q})=c_{1}$, the variable $\dot{q}$ in terms of $t$,
$q$ and the constant of motion $c_{1}$ as follows:
$\dot{q}=\chi_{1}(t,q;c_{1})$ where $\chi_{1}$ is a known function of $t$, $q$
and $c_{1}$. Then, by replacing this function in the second equation
$I_{2}(t,q,\dot{q})=c_{2}$, one obtains the following equation:

$c_{2}=I_{2}\big{(}t,q,\chi_{1}(t,q;c_{1})\big{)}:=\Sigma(t,q;c_{1}),$ (2.3)

where $\Sigma$ is a known function of variables $t$, $q$ and constant of
motion $c_{1}$. Equation (2.3) can be solved algebraically and obtained
$q=\Pi(t;c_{1},c_{2})$, where $\Pi$ is a known function of variable $t$ and
constants of motion $c_{1}$ and $c_{2}$. In this way, the solution of the
problem is complete, and $q=\Pi(t;c_{1},c_{2})$ derived by the ISG-method is
the general solution of the DE (2.1).

Now, to get into the details of the ISG-method, how can one find two
independent invariants $I_{1}(t,q,\dot{q})=c_{1}$ and
$I_{2}(t,q,\dot{q})=c_{2}$? As mentioned before, if the force function
$\phi(t,q,\dot{q})$ in the DE (2.1) is a fractional function of polynomials
like $P$ and $Q$ of variables $t$, $q$ and $\dot{q}$ with coefficients in the
set of complex numbers as $\phi(t,q,\dot{q})=P/Q$, then, to find these two
invariants, one can use the following theorem.

Theorem 1. (Duarte’s integral formula) Let $\ddot{q}=\phi(t,q,\dot{q})$ be the
DE of a particle with one degree of freedom in the configuration space
$\mathbb{Q}=(q)$, where the force function $\phi(t,q,\dot{q})$ is a fractional
function of polynomials $P$ and $Q$ of variables $t$, $q$ and $\dot{q}$. If
this second-order ODE admits a first integral as $I(t,q,\dot{q})=c$, then

$\displaystyle
I(t,q,\dot{q})=r_{1}+r_{2}-\int\Big{[}R+\frac{\partial}{\partial\dot{q}}(r_{1}+r_{2})\Big{]}d\dot{q},$
(2.4)

where

$\displaystyle r_{1}=\int
R(\phi+S\dot{q})dt,~{}~{}~{}~{}~{}r_{2}=-\int\Big{(}RS+\frac{\partial
r_{1}}{\partial q}\Big{)}dq.$ (2.5)

In equations (2.4) and (2.5), $S$ and $R$, which satisfy in the PS determining
equations:

$\displaystyle D[S]$ $\displaystyle=$
$\displaystyle-\phi_{{}_{q}}+S\phi_{{}_{\dot{q}}}+S^{2},$ (2.6) $\displaystyle
D[R]$ $\displaystyle=$ $\displaystyle-R(S+\phi_{{}_{\dot{q}}}),$ (2.7)
$\displaystyle R_{q}$ $\displaystyle=$ $\displaystyle
R_{{}_{\dot{q}}}S+RS_{{}_{\dot{q}}},$ (2.8)

are called the “null form” and the “integrating factor” of the DE (2.1),
respectively [21, 22, 28]. Here, $\phi_{{}_{q}}:=\frac{\partial\phi}{\partial
q},\phi_{{}_{\dot{q}}}:=\frac{\partial\phi}{\partial\dot{q}}$ and so on.

In order to calculate the functionally independent invariants
$I_{1}(t,q,\dot{q})=c_{1}$ and $I_{2}(t,q,\dot{q})=c_{2}$ associated to the DE
(2.1), we must first obtain the null form $S$ and integrating factor $R$,
using the PS determining equations (2.6)-(2.8). Solving these differential
equations for $S$ and $R$ is difficult, except for certain simple cases.
Therefore, to obtain these two basic functions in the PS procedure, we must
resort to another symmetry methods. These methods are Lie point symmetry,
$\lambda$-symmetry and DPs.

If the calculation of functions $S$ and $R$ is not achieved by the PS
determining equations, then one must do them indirectly. For this purpose,
suppose [27]

$\displaystyle\mathbf{\Phi}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2},$
$\displaystyle(t,q)\mapsto(\overline{t},\overline{q})=\mathbf{\Phi}(t,q;\epsilon)=\big{(}T(t,q;\epsilon),Q(t,q;\epsilon)\big{)},$
(2.9)

with the transformation equations

$\displaystyle{t}\rightarrow{\bar{t}}=t+\epsilon{\tau}(t,q)+{\cal
O}(\epsilon^{2}),$ $\displaystyle
q\rightarrow{\bar{q}}={q}+\epsilon{\xi}(t,{q})+{\cal
O}(\epsilon^{2}),~{}~{}\epsilon\ll 1,$ (2.10)

be the group of one-parameter Lie point transformations such that the DE (2.1)
under which remains invariant. In equations (2), functions $\tau$ and $\xi$
defined by

$\displaystyle{\tau}(t,{q})$ $\displaystyle=$ $\displaystyle\frac{\partial
T(t,{q};\epsilon)}{\partial\epsilon}{|_{{}_{\epsilon=0}}},$ (2.11)
$\displaystyle{\xi}(t,{q})$ $\displaystyle=$
$\displaystyle\frac{\partial{Q}(t,{q};\epsilon)}{\partial\epsilon}{|_{{}_{\epsilon=0}}},$
(2.12)

are called infinitesimals of the group of transformations $\mathbf{\Phi}$. By
these two functions, the infinitesimal generator vector of the group of
transformations $\mathbf{\Phi}$ is defined as

$\displaystyle{\bf X}=\tau(t,q)\frac{\partial}{\partial
t}+\xi(t,q)\frac{\partial}{\partial q}.$ (2.13)

In order to find the Lie point symmetries associated to the DE (2.1), we must
first obtain the infinitesimals $\tau$ and $\xi$. For this purpose, we use the
following theorem.

Theorem 2. (Lie’s invariance condition) Suppose for a dynamical system with
one degree of freedom, the second-order ODE
$H(t,q,\dot{q},\ddot{q}):=\ddot{q}-\phi(t,q,\dot{q})=0$ be the DE (or
governing equation) of the particle in the configuration space, and let ${\bf
X}=\tau(t,x)\partial_{{}_{t}}+\xi(t,q)\partial_{{}_{q}}$ be the infinitesimal
generator of the one-parameter Lie group of transformations (2) with
transformation equations (2) acting on the $(t,q)$-space. The group of
transformations $\mathbf{\Phi}$ is admitted by the DE (2.1), if and only if
the following condition is holds:

$\displaystyle{{{\bf X}{{}^{(2)}}H(t,q,\dot{q},\ddot{q})}{\Big{|}}}_{H=0}=0,$
(2.14)

where

$\displaystyle{\bf X}{{}^{(2)}}:=\tau\frac{\partial}{\partial
t}+\xi\frac{\partial}{\partial
q}+(\dot{\xi}-\dot{q}\dot{\tau})\frac{\partial}{\partial\dot{q}}+(\ddot{\xi}-2\ddot{q}\dot{\tau}-\dot{q}\ddot{\tau})\frac{\partial}{\partial\ddot{q}},$
(2.15)

is the 2th extended (prolonged) infinitesimal generator of the vector field
${\bf X}$. We note that equation (2.14) is called the Lie’s invariance
condition which is the fundamental equation of the Lie symmetry analysis [27,
29, 30].

In the other words, the establishing equation (2.14) is a necessary and
sufficient condition for the ODE $H=0$ to admit $\mathbf{\Phi}$ with
infinitesimal generator $\bf X$ as a one-parameter Lie group of
transformations. By solving the Lie’s invariance condition, which are in fact
a system of second-order partial differential equations for the unknown
functions $\tau$ and $\xi$, all infinitesimal generators $\bf X$ associated to
the ODE $H=0$ can be found.

Assume that the ODE $H=0$ has at least two Lie point symmetries and
${\Phi_{1}}$ and ${\Phi_{2}}$ be arbitrarily two of them with the
infinitesimal generators ${\bf
X}_{1}=\tau_{1}\partial_{{}_{t}}+\xi_{1}\partial_{{}_{q}}$ and ${\bf
X}_{2}=\tau_{2}\partial_{{}_{t}}+\xi_{2}\partial_{{}_{q}}$, respectively. By
the symmetries ${\bf X}_{1}$ and ${\bf X}_{2}$, we define the characteristics
$Q_{1}:=\xi_{1}-\dot{q}\tau_{1}$ and $Q_{2}:=\xi_{2}-\dot{q}\tau_{2}$,
respectively. Similarly, by the characteristics $Q_{1}$ and $Q_{2}$ we define
the functions $\lambda_{i}:=D[Q_{i}]/Q_{i}$, ($i$=1, 2). It can be shown that
$-\lambda_{i}$’s are solutions of the determining equation (2.6). Thus, we
conclude that the null forms are $S_{i}=-\lambda_{i}$. Now, by using the
Darboux’s eigenvalue equation444In this equation, $F$’s are DPs associated to
the DE (2.1), and the eigenvalue
$K(t,q,\dot{q})\big{(}\phi_{\dot{q}}(t,q,\dot{q})\big{)}$ is called the
cofactor of $F$ when DE is an explicit function of $t$ (when DE is not an
explicit function of $t$) and has degree at most 2. [21, 27]:

$\displaystyle D[F]=\begin{cases}\phi_{{}_{\dot{q}}}F&~{}~{}~{}~{}$if DE is
not an explicit function of$~{}t,\\\ KF&~{}~{}~{}~{}$if DE is an explicit
function of$~{}t.\\\ \end{cases}$ (2.16)

one can obtain the DPs $F$ associated to the DE (2.1). It can be shown that
the ratio $Q/F$ is a solution of the determining equation (2.7), and hence
$R=Q/F$. Therefore, two 2-tuples $(S_{1},R_{1})$ and $(S_{2},R_{2})$ are
obtained as follow:

$\displaystyle(S_{1},R_{1})=(-\lambda_{1},\frac{Q_{1}}{F}),~{}~{}~{}~{}~{}(S_{2},R_{2})=(-\lambda_{2},\frac{Q_{2}}{F}).$
(2.17)

Thus, by substituting each of these 2-tuples $(S_{1},R_{1})$ and
$(S_{2},R_{2})$ in the Duarte’s integral formula (2.4) one can obtain the
first integrals corresponding to each of these 2-tuples. Certainly, these
first integrals are independent of each other. Therefore, they are what we
needed to solve the DE (2.1). It should be noted if the DE (2.1) does not have
any Lie point symmetry that can be used to calculate the null form $S$ and the
integrating factor $R$, then, one must use its generalized symmetries. The
most suitable of these symmetries is $\lambda$-symmetry. Similar to Lie point
symmetry, $\lambda$-symmetry has also a symmetry condition. If this symmetry
condition is used for the DE (2.1), then, the function $\lambda(t,q,\dot{q})$
and its infinitesimals $\tau$ and $\xi$ are obtained. When the functions
$\lambda$, $\tau$ and $\xi$ are known, then the first integral can be
calculated by a process which is given in Ref. [27].

##  3 FLRW cosmological model in BD theory

The purpose of this section is to study the FLRW cosmological model in the BD
theory. Before proceeding to do this, let us give a short review of the theory
of BD.

###  3.1 A short review of the theory of BD

The BD theory of gravity is a generalization of Einstein’s theory of general
relativity, which was formulated by Brans and Dicke [3, 4]. As mentioned
earlier, this theory is consistent with Mach’s principle (the stronger form)
and relies less on the inherent properties of space. In this theory, in
addition to the tensor field that describes the geometry of space-time, a
scalar field is also included. Hence, BD theory is also called scalar-tensor
theory. To obtain the field equations of this theory, let us consider
Einstein-Hilbert action in the theory of GR

$\displaystyle I_{{}_{EH}}=\int_{\cal M}\sqrt{-g}\Big{(}-\frac{\cal R}{16\pi
G}+{\mathcal{L}}_{{}_{M}}\Big{)}d^{4}x.$ (3.1)

where ${\cal R}$ is the Ricci scalar associated to the four-dimensional space-
time manifold ${\cal M}$ with the local coordinates
$x^{\mu}=(x^{0},x^{1},x^{2},x^{3})$, and
${\mathcal{L}}_{{}_{M}}={\mathcal{L}}_{{}_{M}}(g_{{}_{\mu\nu}},\partial_{{}_{k}}g_{{}_{\mu\nu}})$
is the Lagrangian density of the ordinary matter and
$g:=det~{}g_{{}_{\mu\nu}}$. This action in the BD theory and in the “Jordan
frame” is generalized as follows [10, 12]:

$\displaystyle I_{{}_{BD}}=\frac{1}{16\pi}\int_{\cal
M}\sqrt{-g}\Big{(}\varphi{\cal
R}-\frac{\omega}{\varphi}~{}g^{\mu\nu}~{}\nabla_{\mu}\varphi~{}\nabla_{\nu}\varphi-V(\varphi)+16\pi~{}{\mathcal{L}}_{{}_{M}}\Big{)}d^{4}x,$
(3.2)

where $\varphi$ is the BD scalar field, $\omega$ is the adjustable BD coupling
parameter, $V(\varphi)$ is the potential energy of the field $\varphi$ and
$\nabla_{\mu}$ denotes the covariant derivative operator in the space-time
$x^{\mu}$. Moreover, the BD scalar field $\varphi$ inversely proportional to
the effective gravitational constant $G_{{}_{eff}}$ by the following relation

$\displaystyle\varphi=\frac{1}{G_{eff}}~{}\frac{4+2\omega}{3+2\omega},$ (3.3)

where $G_{eff}$ is equal to the constant $G$ in the limit when $\omega$ tends
to infinity. In addition, ${\mathcal{L}}_{M}$ as the Lagrangian density is
minimally coupled to the BD scalar field $\varphi$. In analogy to the GR
theory, the law of conservation of energy-momentum of the ordinary matter
$\nabla_{\mu}T^{\mu\nu}_{M}=0$ and the definition of the energy-momentum
tensor (EMT) of the ordinary matter fields in the four-dimensional space-time
is given by

$\displaystyle T^{\mu\nu}_{M}:=\frac{2}{\sqrt{-g}}~{}\frac{\partial}{\partial
g_{{}_{\mu\nu}}}~{}(\sqrt{-g}~{}\mathcal{L}_{{}_{M}}).$ (3.4)

By varying the action (3.2) with respect to the scalar field $\varphi$ and
metric $g^{\mu\nu}$ one can get the following equations, respectively [31]

$\displaystyle\frac{2\omega}{\varphi}~{}\square\varphi+{\cal
R}-\frac{\omega}{\varphi^{2}}~{}\nabla^{\alpha}\varphi~{}\nabla_{\alpha}\varphi-\frac{dV}{d\varphi}=0,$
(3.5) $\displaystyle{\cal R}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}~{}{\cal
R}=\frac{8\pi}{\varphi}~{}T_{M\mu\nu}+\frac{\omega}{\varphi^{2}}(\nabla_{\mu}\varphi~{}\nabla_{\nu}\varphi-\frac{1}{2}g_{\mu\nu}~{}\nabla^{\alpha}\varphi~{}\nabla_{\alpha}\varphi)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{1}{\varphi}(\nabla_{{}_{\mu}}\nabla_{{}_{\nu}}\varphi-
g_{{}_{\mu\nu}}\square\varphi)-\frac{V}{2\varphi}g_{{}_{\mu\nu}},$ (3.6)

where ${\cal R}_{\mu\nu}$ is the Ricci tensor and
$\square:=\nabla^{\alpha}~{}\nabla_{\alpha}$ is the covariant d’Alemberian
operator of the metric. Note that equation (3.5) is known as the “modified
Klein-Gordon equation”, and equation (3.6) is a generalization of the EFEs
${\cal R}_{\mu\nu}-(1/2)g_{\mu\nu}~{}{\cal R}+\Lambda
g_{\mu\nu}=8\pi~{}G~{}T_{M\mu\nu}$. Now, by contracting on equation (3.6) one
obtains

$\displaystyle{\cal R}=-\frac{8\pi
T_{M}}{\varphi}+\frac{\omega}{\varphi^{2}}~{}\nabla^{\alpha}\varphi~{}\nabla_{\alpha}\varphi+\frac{3\square\varphi}{\varphi}+\frac{2V}{\varphi}.$
(3.7)

In order to eliminate the Ricci scalar ${\cal R}$ from (3.5) one may insert
(3.7) into (3.5). It then results

$\displaystyle\square\varphi=\frac{8\pi}{3+2\omega}~{}(8\pi
T_{M}+\varphi\frac{dV}{d\varphi}-2V),$ (3.8)

where $T_{M}:=g^{\mu\nu}~{}{T_{{}_{M}}}_{\mu\nu}$ is the trace of the energy-
momentum tensor of the ordinary matter fields. The equation (3.8) together
with (3.6) form a system of coupled second-order non-linear differential
equations which are called the general form of the BD equations. Choosing the
potential $V(\varphi)$ in the form $V(\varphi)=2\Lambda\varphi$, the field
equations (3.6) and (3.8) yield the following standard BD equations

$\displaystyle\square\varphi$ $\displaystyle=$
$\displaystyle\frac{8\pi}{3+2\omega}~{}T_{M}-\frac{2\Lambda\varphi}{3+2\omega},$
(3.9) $\displaystyle{\cal R}_{\mu\nu}-\frac{1}{2}g_{\mu\nu}~{}{\cal R}+\Lambda
g_{\mu\nu}$ $\displaystyle=$
$\displaystyle\frac{8\pi}{\varphi}~{}T_{M\mu\nu}+\frac{\omega}{\varphi^{2}}(\nabla_{\mu}\varphi~{}\nabla_{\nu}\varphi-\frac{1}{2}g_{\mu\nu}~{}\nabla^{\alpha}\varphi~{}\nabla_{\alpha}\varphi)$
(3.10)
$\displaystyle+\frac{1}{\varphi}(\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\square\varphi).$

This system of the field equations, without considering the cosmological
constant $\Lambda$, were obtained for the first time by Brans and Dicke in
1961, and known as the BD equations in the presence of the cosmological
constant $\Lambda$. These equations form the basis of the BD theory. In fact,
BD theory with the BD equations as the field equations is an extended gravity
theory of the “scalar-tensor” type, in which the potential energy function
sets to $V(\varphi)=2\Lambda\varphi$. In BD theory, although the initial value
of the BD coupling parameter was very small, its present value exceeds 500 in
many cosmological models, such that $\omega>500$ is in good agreement with
experimental tests and general relativity predictions (for example, neutron
stars, black holes, gravitational waves, binary pulsars, etc.) [32]. The BD
theory of gravity reduces the GR theory in the $\omega\rightarrow\infty$ limit
[33, 34].

In the following, we will use the ISG-method to solve the DEs of the FLRW
cosmological model in the framework of the BD theory. In fact, we will obtain
the cosmological solutions of the BD equations (3.9) and (3.10) by the ISG-
method.

###  3.2 FLRW cosmological model in BD theory

In this subsection, the FLRW cosmological model is briefly studied in the
framework of the BD theory. We consider a homogeneous and isotropic universe
as a cosmological model for our study. To a good approximation, this universe
is described by the FLRW metric. This metric in the coordinates
$x^{\mu}=(t,r,\theta,\varphi)$ is defined as

$\displaystyle
ds^{2}=-dt^{2}+a^{2}(t)\Big{(}\frac{dr^{2}}{1-kr^{2}}+r^{2}~{}d\theta^{2}+r^{2}~{}sin^{2}\theta~{}d\varphi^{2}\Big{)},$
(3.11)

where $a(t)$ is the scale factor, and the number $k$ can admits three values
of $-1,1,0$ for spaces of negative, positive and zero curvature, which
represents closed, open and flat universes, respectively [35, 36]. The Weyl’s
principle requires that the constituent material of a homogeneous and
isotropic universe is a perfect fluid on large scales with high accuracy. So,
according to the Weyl’s principle, the FLRW universe can be considered as a
perfect fluid with a good approximation. If the energy density of this fluid
be $\rho c^{2}$ and its pressure $p$, then its normal barotropic equation of
state is

$\displaystyle p=w\rho c^{2},$ (3.12)

where $w$ is called the perfect fluid state parameter and $c$ is the velocity
of light in vacuum. The equation of state (3.12) includes most of the
interesting states that are important in cosmology. For example, for dust:
$w=0$, for radiation: $w=1/3$, for false vacuum: $w=-1$, and for stiff fluid:
$w=1$ [37, 38]. In the present work, a system of units is used in which the
velocity of light in vacuum is set to unity. Furthermore, the components of
the EMT of the fluid will be as follows:

$\displaystyle T_{M\mu\nu}=(\rho+p)U_{\mu}~{}U_{\nu}+pg_{\mu\nu},$ (3.13)

where $U^{\mu}=(1,0,0,0)$ is the co-moving 4-velocity. The trace of this
tensor is given by

$\displaystyle T_{M}=3p-\rho.$ (3.14)

The non-vanishing components of the BDE (3.10) read

$\displaystyle\frac{3\dot{a}^{2}}{a^{2}}-\Lambda$ $\displaystyle=$
$\displaystyle\frac{8\pi}{\varphi}\rho+\frac{\omega}{2}\frac{\dot{\varphi}^{2}}{\varphi^{2}}-3\frac{\dot{a}}{a}\frac{\dot{\varphi}}{\varphi},$
(3.15) $\displaystyle-\frac{2\ddot{a}}{a}-\frac{\dot{a}^{2}}{a^{2}}+\Lambda$
$\displaystyle=$
$\displaystyle\frac{8\pi}{\varphi}p+\frac{\omega}{2}\frac{\dot{\varphi}^{2}}{\varphi^{2}}+2\frac{\dot{a}}{a}\frac{\dot{\varphi}}{\varphi}+\frac{\ddot{\varphi}}{\varphi}.$
(3.16)

Also, using equation (3.14), one concludes that equation (3.9) becomes

$\displaystyle\frac{\ddot{\varphi}}{\varphi}+3\frac{\dot{a}}{a}\frac{\dot{\varphi}}{\varphi}=\frac{2\Lambda}{3+2\omega}+\frac{8\pi}{\varphi}\Big{(}\frac{\rho-3p}{3+2\omega}\Big{)}.$
(3.17)

Note that equation (3.17) is not an independent equation in BD theory, because
it can be obtained by using Bianchi identities $\nabla_{\mu}~{}G^{\mu\nu}=0$
and equations (3.15) and (3.16). So, we have only two independent equations to
determine four unknown functions $a(t)$, $\varphi(t)$, $\rho(t)$ and $p(t)$.
As such this system of equations does not seem to have a unique solution.
Hence, to find these unknowns uniquely, two equations are needed. One of these
is the equation of state of the universe, (3.12). To obtain the second
equation we use the following Corollary.

Corollary1. In the spatially flat $(k=0)$ FLRW cosmological model in the
framework of the BD theory, the necessary and sufficient condition for the
deceleration parameter of the universe $q:=-a\ddot{a}/\dot{a}^{2}$ to be a
constant, is that a power-law between the cosmic scale factor $a(t)$ and the
Brans-Dicke scalar field $\varphi(t)$ holds in the following form [35, 36, 37,
38, 39, 40]:

$\displaystyle\varphi a^{n}=C,$ (3.18)

where $n$ and $C$ are the functions of the parameters $w$, $\omega$ and
constants $G$ and $\Lambda$, that is, $C:=C(w,\omega,G,\Lambda)$ and
$n:=n(w,\omega,G,\Lambda)$.
Proof. See Ref. [39].

The history of this Corollary backs to the Dirac’s hypothesis that states the
constant $G$ is a variable quantity and should be related by a power-law
relation such as $G=\acute{C}a^{n}$ to the scale factor of the universe, where
$\acute{C}$ and $n$ are some constants. Since in the BD theory, the scalar
field $\varphi$ is proportional to the inverse of the Newton’s gravitational
constant, that is, $\varphi=c_{{}_{0}}/G$, (where $c_{{}_{0}}$ is the
proportionality constant), combining this relation with $G={C^{\prime}}a^{n}$
yields the power-law in the form (3.18) in which $C:=c_{{}_{0}}/{C^{\prime}}$.
By combining equations (3.15)-(3.17) with equation of state (3.12) and power-
law relation (3.18), we then obtain

$\displaystyle\frac{8\pi\rho}{C}a^{n}$ $\displaystyle=$
$\displaystyle\frac{\dot{a}^{2}}{a^{2}}\Big{(}\frac{6-\omega
n^{2}-6n}{2}\Big{)}-\Lambda,$ (3.19) $\displaystyle\frac{8\pi w\rho}{C}a^{n}$
$\displaystyle=$
$\displaystyle(n-2)\frac{\ddot{a}}{a}+\frac{\dot{a}^{2}}{a^{2}}\Big{[}\frac{2(n-1)-(\omega+2)n^{2}}{2}\Big{]}+\Lambda,$
(3.20) $\displaystyle\frac{8\pi}{C}a^{n}(\rho-3w\rho)$ $\displaystyle=$
$\displaystyle(3+2\omega)\Big{[}-n\frac{\ddot{a}}{a}+n(n-2)\frac{\dot{a}^{2}}{a^{2}}\Big{]}-2\Lambda.$
(3.21)

The combination of these three differential equations after performing some
algebraic calculations leads to the following non-linear second-order ODE for
the scale factor $a(t)$

$\displaystyle 2(\omega n+3)\frac{\ddot{a}}{a}+(6+4n\omega-\omega
n^{2})\frac{\dot{a}^{2}}{a^{2}}=2\Lambda,~{}~{}~{}~{}~{}~{}~{}~{}~{}\omega
n+3\neq 0.$ (3.22)

This equation can be rewritten in the following form

$\displaystyle\ddot{a}=\alpha\frac{\dot{a}^{2}}{a}+\beta a,$ (3.23)

where the new parameters $\alpha$ and $\beta$ in terms of the old parameters
$\omega$, $w$, $n$ and the constant $\Lambda$ are defined by

$\displaystyle\alpha:=-\frac{6+4n\omega-\omega n^{2}}{2(\omega
n+3)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\beta:=\frac{\Lambda}{\omega
n+3}.$ (3.24)

Therefore, the combination of BD equations with the fluid state equation of
the universe and the power-law relation leads to a second-order ODE for the
cosmic scale factor $a(t)$ in the form (3.23). In this work, we call this
equation, the “Friedman-Brans-Dicke equation” (FBD-equation). The FBD-equation
alone describes the dynamics of spatially flat $(k=0)$ FLRW universe in the
framework of the BD theory. Equation (3.23) can be imagined as the equation of
motion of a particle with a unit mass in one-dimension. By defining the cosmic
scale factor as the coordinate of this particle, i.e., $q(t):=a(t)$, equation
(3.23) is the familiar form of Newton’s second law

$\displaystyle\ddot{q}=\phi(t,q,\dot{q}),$ (3.25)

where

$\displaystyle\phi(t,q,\dot{q}):=\alpha\frac{\dot{q}^{2}}{q}+\beta q,$ (3.26)

is the component of the force function acting on the particle, and $\ddot{q}$
is the component of its acceleration vector. Thus, to analyze the dynamics of
cosmological model, it is enough to solve only the FBD-equation by the ISG-
method. We will do this in the next section.

##  4 Solving the FBD-equation

In this section we employ the ISG-method to solve the FBD-equation (3.25). To
this end, we use the following steps:

(1)-We re-consider spatially flat $(k=0)$ FLRW cosmological model in the
framework of BD theory with a dynamical system $S^{2}_{1}$. As mentioned
earlier, the governing equation of this system is
$\ddot{q}=\phi(t,q,\dot{q})$, where $\phi(t,q,\dot{q})$ is the force function
of a particle that moves in a $1$-dimensional mini-super space with
configuration $\mathbb{Q}=(q)$ where $q:=a(t)$ is the coordinate of the
particle.

(2)-A one-parameter Lie group of point transformations
${\bf\Phi}(t,q;\varepsilon)$ with the transformation equations (2) must be
found (if there is a solution) such that the FBD-equation (3.25) remains
invariant under this group of transformations. For this purpose, by using
Lie’s invariance condition, a set of the PDEs should be extracted for the
infinitesimals $\tau$ and $\xi$ of group ${\bf\Phi}$. By simultaneously
solving this system of PDEs, the infinitesimals $\tau$ and $\xi$ as functions
of $t$ and $q$ should be obtained: $\tau=\tau(t,q)$, $\xi=\xi(t,q)$. For our
case, these functions are obtained as $\tau=\tau(t,q)=c_{1}$,
$\xi=\xi(t,q)=c_{2}q$ where $c_{1}$ and $c_{1}$ are some arbitrary constants
[27]. Then, by using these results one must find all Lie point symmetry
vectors of the transformations group ${\bf\Phi}$ such that for our case these
vectors are obtained to be $\mathbf{X_{1}}=\partial_{t}$ and
$\mathbf{X_{2}}=q\partial_{q}$ which are respectively the infinitesimal
generator of the translation group along the time
${\bf\Phi}_{1}(t,q;\varepsilon)=(t+\varepsilon,q)$ and the scaling group
${\bf\Phi}_{2}(t,q;\varepsilon)=(t,e^{\varepsilon}q)$ [27].

(3)-To calculate the first integral corresponding to the DE
$\ddot{q}=\phi(t,q,\dot{q})$ of the particle in the configuration space
$\mathbb{Q}=(q)$, which is equation (2.4) of Theorem 1, the null form $S$ and
integrating factor $R$ should be obtained by using the determining equations
(2.6)-(2.8). It is difficult to solve these equations simultaneously, except
in special simple cases. Therefore, to obtain these two basic functions in the
extended PS method, we must resort to other symmetry methods. The most
important other symmetry methods, that can be used here, are: $(a)$ Lie point
symmetry, $(b)$ $\lambda$-symmetry, $(c)$ DPs. Therefore, the functions $S$
and $R$ should be calculated indirectly by the symmetry methods $(a)$-$(c)$.

(4)-Using the Lie point symmetry mentioned in section 2, the characteristic
$Q(t,q,\dot{q}):=\xi-\dot{q}\tau$ and then $\lambda(t,q,\dot{q}):=D[Q]/Q$ must
be defined. It can be shown that $-D[Q]/Q$ is a solution of the determining
equation (2.6), in such a way that $S(t,q,\dot{q})=-\lambda(t,q,\dot{q})$.

(5)-By using Darboux’s eigenvalue equation (2.16) [21], the DP
${F}(t,q,\dot{q})$ and the eigenvalue
$K(t,q,\dot{q})(\phi_{\dot{q}}(t,q,\dot{q}))$ corresponding to this polynomial
should be obtained. It can be shown that the ratio $Q/F$ is a general solution
of the determining equation (2.7). So, $R(t,q,\dot{q})=Q/F$. In this way, the
null form $S$ and the integrating factor $R$ are obtained indirectly by the
Lie point symmetry method and DPs, respectively, without the need to solve
their determining equations.

Let us turn into the main goal of this section which is nothing but
calculating the first integrals $I_{1}$ and $I_{2}$ corresponding to the FBD-
equation (3.25). Since the FBD-equation admits two Lie point symmetries, both
$\lambda$-symmetries associated to FBD-equation can be only calculated by
using the relation $\lambda=D[Q]/Q$ without solving the $\lambda$-symmetry
condition [27]. Therefore, both $\lambda$-symmetries associated to the Lie
symmetry vectors $\mathbf{X}_{1}=\partial_{t}$ and
$\mathbf{X}_{2}=q\partial_{q}$ are, respectively,

$\displaystyle\lambda_{1}(t,q,\dot{q})$ $\displaystyle=$
$\displaystyle\frac{D[Q_{1}]}{Q_{1}}=\frac{D[-\dot{q}]}{-\dot{q}}=\alpha\frac{\dot{q}}{q}+\beta\frac{q}{\dot{q}},$
(4.1) $\displaystyle\lambda_{2}(t,q,\dot{q})$ $\displaystyle=$
$\displaystyle\frac{D[Q_{2}]}{Q_{2}}=\frac{D[q]}{q}=\frac{\dot{q}}{q}.$ (4.2)

It can be easily check that functions $-D[Q_{i}]/{Q_{i}}$ are solutions of the
determining equations (2.6). Thus, for the FBD-equation, these functions are
nothing but two null forms associated to the Lie symmetry vectors
$\textbf{X}_{1}=\partial_{t}$ and $\textbf{X}_{2}=q\partial_{q}$:

$\displaystyle S_{1}(t,q,\dot{q})$ $\displaystyle=$
$\displaystyle-\lambda_{1}(t,q,\dot{q})=-\alpha\frac{\dot{q}}{q}-\beta\frac{q}{\dot{q}},$
(4.3) $\displaystyle S_{2}(t,q,\dot{q})$ $\displaystyle=$
$\displaystyle-\lambda_{2}(t,q,\dot{q})=-\frac{\dot{q}}{q}.$ (4.4)

On the other hand, the determining equation for DPs of the FBD-equation
$\ddot{q}=\phi(t,q,\dot{q})$ is (2.16). As mentioned above, $Q/F$ is a
solution of the determining equation (2.7) such that $R=Q/F$, and hence by
considering $S=-D[Q]/Q$ we arrive at

$\displaystyle
D[\frac{Q}{F}]+\frac{Q}{F}\Big{(}-\frac{D[Q]}{Q}+\phi_{\dot{q}}\Big{)}=0.$
(4.5)

Furthermore, the characteristics $Q=\xi-\dot{q}\tau$ obtained from the
$\lambda$-symmetry and the DPs $F$ obtained from solving of the Darboux
eigenvalue equation associated to the FBD-equation gives us the integrating
factor $R$ as $R=Q/F$. Therefore, by using $\lambda$-symmetry and DPs without
solving the determining equations, we can obtain $2$-tuple $(S,R)$ as follows
$(S,R)=(-D[Q]/Q,Q/F)$. In the ISG-method, the basic quantities of the Lie
point symmetry, extended PS method, $\lambda$-symmetry, and DPs read $\tau$,
$\xi$, $\lambda$, $F$, $R$ and $S$. These quantities are related to each other
as shown in [27].

Now, we consider one of the Lie symmetries, for example
$\textbf{X}_{1}=\partial_{t}$. This vector is the infinitesimal generator of
the group of translation time:
${\bf\Phi}_{1}(t,q;\varepsilon)=(t+\varepsilon,q)$. The characteristic
corresponding to this symmetry vector is
$Q_{1}=\xi_{1}-\dot{q}\tau_{1}=-\dot{q}$ and $\lambda$ associated to this
characteristic is given by (4.1). Therefore, the null form associated to the
infinitesimal generator $\textbf{X}_{1}=\partial_{t}$ is obtained from
equation (4.3). We find the integrating factor $R_{1}$ corresponding to the
null form $S_{1}$. For this purpose, by using the characteristic
$Q_{1}=-\dot{q}$ and by considering DP555According to Table 1 of Ref. [27], by
the characteristics $Q_{1}(t,q,\dot{q})=-\dot{q},Q_{2}(t,q,\dot{q})=q$ and DPs
$F_{1}^{(D)},F_{2}^{(D)},\cdots,F_{6}^{(D)}$, two null forms
$S_{i}:=-D[Q_{i}]/Q_{i},(i=1,2)$ and twelve integrating factors
$R_{ij}:=Q{i}/F_{j}^{(D)},(i=1,2,j=1,2,...,6)$ can be defined.
$F_{6}^{(D)}=q^{2\alpha}$, the integrating factor $R_{1}$ reads

$\displaystyle R_{1}=\frac{Q_{1}}{F_{6}^{(D)}}=-\frac{\dot{q}}{q^{2\alpha}},$
(4.6)

and $2$-tuple $(S_{1},R_{1})$ is then obtained to be of the form

$\displaystyle(S_{1},R_{1})=(-\alpha\frac{\dot{q}}{q}-\beta\frac{q}{\dot{q}}~{},~{}-\frac{\dot{q}}{q^{2\alpha}}).$
(4.7)

By substituting this $2$-tuple in the Duarte’s integral formula (2.4), and
then by calculating the integrals, one can get the first integral associated
to the infinitesimal generator $\mathbf{X_{1}}=\partial_{t}$, giving us

$\displaystyle I_{1}(q,\dot{q})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\dot{q}^{2}~{}q^{-2\alpha}-\frac{\beta}{2(1-\alpha)}q^{-2\alpha+2}=c_{1},$
(4.8)

where $c_{1}$ is an arbitrary constant. In the same way, the $2$-tuple
$(S_{2},R_{2})$ corresponding to the symmetry vector
$\textbf{X}_{2}=q\partial_{q}$ can be obtained as follows:

$\displaystyle(S_{2},R_{2})$ $\displaystyle=$
$\displaystyle(-\lambda_{2},Q_{2}/F_{1}^{(D)})=\big{(}-\frac{\dot{q}}{q}~{},~{}\frac{q}{\frac{\beta}{\alpha-1}q^{2}+\dot{q}^{2}}\big{)}.$
(4.9) $\displaystyle I_{2}(t,q,\dot{q})$ $\displaystyle=$
$\displaystyle\begin{cases}(\alpha-1)t+\sqrt{\frac{\alpha-1}{\beta}}\tan^{-1}~{}(\frac{q}{\dot{q}}\sqrt{\frac{\beta}{\alpha-1}})=c_{2},~{}&~{}~{}\frac{\beta}{\alpha-1}>0,{}{}{}{}{}\\\
(\alpha-1)t+\frac{q}{\dot{q}}=c_{2},~{}&~{}~{}\frac{\beta}{\alpha-1}=0,\\\
(\alpha-1)t+\sqrt{|\frac{\alpha-1}{\beta}|}\tanh^{-1}~{}(\frac{q}{\dot{q}}\sqrt{|\frac{\beta}{\alpha-1}|})=c_{2},~{}&~{}~{}\frac{\beta}{\alpha-1}<0.\end{cases}$
(4.10)

According to the above results and as earlier mentioned, a set of twelve
members of the null forms and integrating factors can be constructed, which we
call the PS set:

$\displaystyle
S_{PS}=\\{(S_{i},R_{ij})\\}_{i=1,2,j=1,...,6}=\\{(S_{1},R_{11}),(S_{1},R_{12}),(S_{1},R_{13}),(S_{1},R_{14}),(S_{1},R_{15}),$
$\displaystyle(S_{1},R_{16}),(S_{2},R_{21}),(S_{2},R_{22}),(S_{2},R_{23}),(S_{2},R_{24}),(S_{2},R_{25}),(S_{2},R_{26})\\},$
(4.11)

By using the formula (2.4), a first integral can be constructed with each of
the $2$-tuple of this set. Therefore, for twelve first integrals we have

$\displaystyle
I_{i,ij}=r_{1;ij}+r_{2;ij}-\int\big{[}R_{ij}+\frac{\partial}{\partial\dot{q}}(r_{1;ij}+r_{2;ij})\big{]}d\dot{q}=c_{i,ij},~{}~{}~{}~{}~{}~{}i=1,2;j=1,...,6,$
(4.12)

where $c_{i,ij}$’s are some constants on the solutions of the DE
$\ddot{q}=\phi(t,q,\dot{q})$. Also, according to (2.5), $r_{1;ij}$’s and
$r_{2;ij}$’s are defined as

$\displaystyle r_{1;ij}(t,q,\dot{q})$ $\displaystyle=$ $\displaystyle\int
R_{ij}(\phi+\dot{q}S_{i})dt,$ (4.13) $\displaystyle r_{2;ij}(t,q,\dot{q})$
$\displaystyle=$
$\displaystyle-\int[R_{ij}S_{i}+\frac{\partial}{\partial_{q}}(r_{i;ij})]dq.$
(4.14)

The FBD-equation $\ddot{q}=\phi(t,q,\dot{q})$ admits all members of the PS set
(4.10) as $2$-tuple $(S,R)$. The members of $S_{PS}$ and their corresponding
first integrals have been listed in Table 2 in Ref. [27]. Among the first
integrals of Table 2, only two of them, $I_{1,16}$ and $I_{2,21}$, which we
denoted by $I_{1}$ and $I_{2}$, respectively, are independent and the rest are
dependent on two first integrals $I_{1}$ and $I_{2}$. In the PS set (4.11),
$(S_{1},R_{16})$ and $(S_{2},R_{21})$ are actually the $2$-tuples
$(S_{1},R_{1})$, $(S_{2},R_{2})$, respectively, and the first integrals
associated to these $2$-tuples are $I_{1}=I_{1,16}$ and $I_{2}=I_{2,21}$,
respectively.

As was mentioned earlier, the infinitesimal generator of the translation group
is the symmetry vector $\textbf{X}_{1}=\partial_{t}$. Therefore, the first
integral associated to this symmetry vector gives the energy of the dynamical
system $S_{1}^{2}$. Thus, the invariant $I_{1}$ is the energy of the dynamical
system. Now, consider the independent first integrals (4.8) and (4.10). We
look at these two independent invariants as a system of algebraic equations.
To solve this system of equations we must remove the variable $\dot{q}$ among
the equations of this system. To do this, we find $\dot{q}$ from the equation
$I_{2}(t,q,\dot{q})=c_{2}$, and then employ $I_{1}(q,\dot{q})=c_{1}$. Finally,
the general solutions of the FBD-equation, which include all three cases
$c<0$, $c=0$ and $c>0$, are worked out

$\displaystyle
q(t)=\begin{cases}\Big{(}\frac{2c_{1}}{|c|}\Big{)}^{\frac{1}{2(1-\alpha)}}~{}\sinh^{\frac{1}{1-\alpha}}\Big{[}(\alpha-1)\sqrt{|c|}(t+\acute{c}_{2})\Big{]},~{}~{}~{}~{}~{}~{}&~{}~{}~{}~{}~{}~{}c<0,{}{}{}{}{}\\\
\Big{(}2c_{1}(\alpha-1)^{2}\Big{)}^{\frac{1}{2(1-\alpha)}}\Big{(}t+{c^{\prime}}_{2}\Big{)}^{\frac{1}{1-\alpha}},~{}~{}~{}~{}~{}~{}&~{}~{}~{}~{}~{}~{}c=0,\\\
\Big{(}\frac{2c_{1}}{c}\Big{)}^{\frac{1}{2(1-\alpha)}}~{}\sin^{\frac{1}{1-\alpha}}\Big{[}(\alpha-1)\sqrt{c}(t+\acute{c}_{2})\Big{]},~{}~{}~{}~{}~{}~{}&~{}~{}~{}~{}~{}~{}c>0.\end{cases}$
(4.15)

where $\alpha$ and $\beta$ are given by (3.24). Moreover,

$\displaystyle
c=\frac{\beta}{\alpha-1},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{c^{\prime}}_{2}$
$\displaystyle:=$ $\displaystyle\frac{c_{2}}{1-{\alpha}}.$ (4.16)

By redefinition ${c^{\prime}}_{1}:=(2c_{1})^{\frac{1}{2(1-\alpha)}}$, the
general solutions (4.15) can be written as follows:

$\displaystyle
q(t)=\begin{cases}c^{\prime}_{1}\Big{(}{\dfrac{-2\Lambda}{s}}\Big{)}^{{}^{-\frac{wn+3}{s}}}\sin^{{}^{\frac{2(wn+3)}{s}}}\Big{[}\dfrac{1}{wn+3}\sqrt{\frac{-\Lambda
s}{2}}(t+c^{\prime}_{2})\Big{]},~{}~{}~{}~{}~{}~{}&~{}~{}~{}~{}~{}~{}\Lambda<0,{}{}{}{}{}\\\
c^{\prime}_{1}\Big{(}\dfrac{s}{2(wn+3)}\Big{)}^{{}^{{\frac{2(wn+3)}{s}}}}\Big{(}t+c^{\prime}_{2}\Big{)}^{{}^{\frac{2(wn+3)}{s}}},~{}~{}~{}~{}~{}~{}&~{}~{}~{}~{}~{}~{}\Lambda=0,\\\
c^{\prime}_{1}\Big{(}{\dfrac{2\Lambda}{s}}\Big{)}^{{}^{-\frac{wn+3}{s}}}\sinh^{{}^{\frac{2(wn+3)}{s}}}\Big{[}\dfrac{1}{wn+3}\sqrt{\frac{\Lambda
s}{2}}(t+c^{\prime}_{2})\Big{]},~{}~{}~{}~{}~{}~{}&~{}~{}~{}~{}~{}~{}\Lambda>0.\end{cases}$
(4.17)

where $s:=12+6n\omega-\omega n^{2}$. As we expected, since the FBD-equation is
an ODE of the second-order, its general solution contains two integration
constants ${c^{\prime}}_{1}$ and ${c^{\prime}}_{2}$. These constants can be
determined by using the initial conditions of the problem, i.e.,
$q(t=0)=q_{0}$ and $\dot{q}(t=0)=\dot{q}_{0}$.

##  5  Completing the solutions

Now, we show that the parameter $n$ in the power-law equation (3.18) with the
parameter of state $w$ in equation (3.12), are consistent, provided that there
is a relationship between them. For this purpose, subtracting both sides of
equations (3.19) and (3.20) we obtain

$\displaystyle\frac{8\pi
a^{n}}{C}\rho(1-w)=\frac{\dot{a}^{2}}{a^{2}}\Big{(}\frac{6-\omega
n^{2}-6n}{2}-\frac{2(n-1)-(\omega+2)n^{2}}{2}\Big{)}-(n-2)\frac{\ddot{a}}{a}-2\Lambda.$
(5.1)

By taking two consecutive times of the total derivative with respect to cosmic
time $t$ from equation (3.18) one obtains

$\displaystyle\frac{\dot{a}}{a}$ $\displaystyle=$
$\displaystyle-\frac{1}{n}\frac{\dot{\varphi}}{\varphi},$ (5.2)
$\displaystyle\frac{\ddot{a}}{a}$ $\displaystyle=$
$\displaystyle\frac{n+1}{n^{2}}\frac{\dot{\varphi}^{2}}{\varphi^{2}}-\frac{1}{n}\frac{\ddot{\varphi}}{\varphi}.$
(5.3)

By substituting (3.18), (5.2) and (5.3) into equation (5.1) and after a little
simplification we arrive at

$\displaystyle\ddot{\varphi}-\frac{3}{n}\frac{\dot{\varphi}^{2}}{\varphi}=\frac{2n}{n-2}\Lambda\varphi+8\pi\rho(1-w)\frac{n}{n-2}.$
(5.4)

In addition, one may insert equations (3.12), (3.18) and (5.2) into (3.17) to
get

$\displaystyle\ddot{\varphi}-\frac{3}{n}\frac{\dot{\varphi}^{2}}{\varphi}=\frac{2\Lambda\varphi}{3+2\omega}+\frac{8\pi}{3+2\omega}\rho(1-3w).$
(5.5)

Comparing equations (5.5) and (5.4) yields the following equation:

$\displaystyle
n=\begin{cases}\frac{1-3w}{\omega(w-1)-1},~{}~{}~{}~{}~{}~{}~{}\Lambda=0,~{}w\in[-1,1],\\\
\\\ -\frac{1}{\omega+1},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Lambda\neq
0,~{}w=0.\end{cases}$ (5.6)

It can be seen that the parameter $n$ in the power-law equation (3.18) cannot
admit any value, but only those values satisfying in (5.6). Thus, one may use
(5.6) to write the general solution (4.17) in terms of the parameters
$\omega$, $w$ and $\Lambda$, and hence by changing $q(t)\rightarrow a(t)$ we
get

$\displaystyle
a(t)=\begin{cases}{c^{\prime}}_{1}\Big{(}-\frac{2\Lambda(\omega+1)^{2}}{(4+3\omega)(3+2\omega)}\Big{)}^{-\frac{\omega+1}{3\omega+4}}\sin^{\frac{2(\omega+1)}{3\omega+4}}\Big{[}\sqrt{-(\frac{3\omega+4}{2\omega+3})\frac{\Lambda}{2}}(t+{c^{\prime}}_{2})\Big{]},~{}~{}~{}\Lambda<0,\\\
\\\
{c^{\prime}}_{1}\Big{(}\frac{4+3\omega(1-w^{2})}{2[\omega(1-w)+1]}\Big{)}^{\frac{2[\omega(1-w)+1]}{4+3\omega(1-w^{2})}}\big{(}t+{c^{\prime}}_{2}\big{)}^{\frac{2[\omega(1-w)+1]}{4+3\omega(1-w^{2})}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Lambda=0,\\\
\\\
{c^{\prime}}_{1}\Big{(}\frac{2\Lambda(\omega+1)^{2}}{(4+3\omega)(3+2\omega)}\Big{)}^{-\frac{\omega+1}{3\omega+4}}\sinh^{\frac{2(\omega+1)}{3\omega+4}}\Big{[}\sqrt{(\frac{3\omega+4}{2\omega+3})\frac{\Lambda}{2}}(t+{c^{\prime}}_{2})\Big{]},~{}~{}~{}~{}~{}~{}~{}~{}\Lambda>0.\end{cases}$
(5.7)

According to (5.6), the above equation is actually the solution of the BD
equations for the scale factor of the dust-dominated universe in both cases
$\Lambda<0$ and $\Lambda>0$, while the solution of the case $\Lambda=0$ can be
included the universe fulfilling with the dust, radiation and false vacuum, as
well as stiff fluid. Having the cosmic scale factor $a(t)$, one can calculate
the BD scalar field $\varphi(t)$. Using equations (3.18), (5.6) and (5.7) and
simplifying them we then get

$\displaystyle\varphi(t)=\begin{cases}C{{c^{\prime}}_{1}}^{\frac{1}{\omega+1}}\Big{(}-\frac{2\Lambda(\omega+1)^{2}}{(4+3\omega)(3+2\omega)}\Big{)}^{-\frac{1}{3\omega+4}}\sin^{\frac{2}{3\omega+4}}\Big{[}\sqrt{-(\frac{3\omega+4}{2\omega+3})\frac{\Lambda}{2}}(t+{c^{\prime}}_{2})\Big{]},~{}~{}~{}~{}\Lambda<0,\\\
\\\
C{{c^{\prime}}_{1}}^{\frac{3w-1}{\omega(w-1)-1}}\Big{(}\frac{4+3\omega(1-w^{2})}{2[\omega(1-w)+1]}\Big{)}^{\frac{2(1-3w)}{4+3\omega(1-w^{2})}}\big{(}t+{c^{\prime}}_{2}\big{)}^{\frac{2(1-3w)}{4+3\omega(1-w^{2})}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Lambda=0,\\\
\\\
C{{c^{\prime}}_{1}}^{\frac{1}{\omega+1}}\Big{(}\frac{2\Lambda(\omega+1)^{2}}{(4+3\omega)(3+2\omega)}\Big{)}^{-\frac{1}{3\omega+4}}\sinh^{\frac{2}{3\omega+4}}\Big{[}\sqrt{(\frac{3\omega+4}{2\omega+3})\frac{\Lambda}{2}}(t+{c^{\prime}}_{2})\Big{]},~{}~{}~{}~{}~{}~{}~{}~{}\Lambda>0.\end{cases}$
(5.8)

In order to calculate the energy density of the universe with the equation of
state $p=\omega\rho$, in all three cases $\Lambda<0$, $\Lambda=0$, and
$\Lambda>0$ we use equations (3.19), (5.6) and (5.7). By doing some algebraic
calculations, we get the following solution for the energy density of the
universe

$\displaystyle\rho(t)=\begin{cases}(-\frac{\Lambda
C}{8\pi}){{c^{\prime}}_{1}}^{-\frac{1}{\omega+1}}\Big{(}-\frac{2\Lambda(\omega+1)^{2}}{(4+3\omega)(3+2\omega)}\Big{)}^{-\frac{1}{3\omega+4}}\sin^{-\frac{6(\omega+1)}{3\omega+4}}\Big{[}\sqrt{-(\frac{3\omega+4}{2\omega+3})\frac{\Lambda}{2}}(t+\acute{c}_{2})\Big{]},~{}~{}\Lambda<0,\\\
\\\
M(\omega,w)X(\omega,w)(t+{c^{\prime}}_{2})^{N(\omega,w)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Lambda=0,\\\
\\\ (\frac{\Lambda
C}{8\pi}){{c^{\prime}}_{1}}^{-\frac{1}{\omega+1}}\Big{(}\frac{2\Lambda(\omega+1)^{2}}{(4+3\omega)(3+2\omega)}\Big{)}^{-\frac{1}{3\omega+4}}\sinh^{-\frac{6(\omega+1)}{3\omega+4}}\Big{[}\sqrt{(\frac{3\omega+4}{2\omega+3})\frac{\Lambda}{2}}(t+\acute{c}_{2})\Big{]},~{}~{}~{}~{}~{}~{}~{}~{}~{}\Lambda>0,\end{cases}$
(5.9)

where $M$, $X$, and $N$ are the functions of the parameters $\omega$ and $w$
which are defined as

$\displaystyle M(\omega,w)$ $\displaystyle:=$
$\displaystyle\frac{C~{}{{c^{\prime}}_{1}}^{\frac{3w-1}{\omega(w-1)-1}}}{4\pi}\Big{(}\frac{4+3\omega(1-w^{2})}{2[\omega(1-w)+1]}\Big{)}^{\frac{2(1-3w)}{4+3\omega(1-w^{2})}},$
(5.10) $\displaystyle X(\omega,w)$ $\displaystyle:=$
$\displaystyle\frac{6(w-1)^{2}~{}\omega^{2}-\omega[12(w-1)+(1-3w)(3w-5)]+12-18w}{[4+3\omega(1-w^{2})]^{2}},$
(5.11) $\displaystyle N(\omega,w)$ $\displaystyle:=$
$\displaystyle\frac{2(1-3w)-8-6(1-w^{2})\omega}{4+3\omega(1-w^{2})}.$ (5.12)

As mentioned earlier, the cases $\Lambda<0$ and $\Lambda>0$ of solution (5.9)
give the energy density for dust dominated universe only, while the case
$\Lambda=0$ of this solution gives the energy density of the perfect fluid-
filled universe with equation of state (3.12). Accordingly and also using the
equation of state $p=w\rho$ and energy density (5.9), the pressure of the
universe is worked out to be

$\displaystyle
p(t)=\begin{cases}0~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Lambda<0,\\\
\\\
wM(\omega,w)X(\omega,w)(t+{c^{\prime}}_{2})^{N(\omega,w)},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Lambda=0,~{}~{}w\in[-1,1]\\\
\\\
0~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\Lambda>0.\end{cases}$
(5.13)

According to this solution, the pressure of the dust universe for $\Lambda<0$,
$\Lambda=0$ and $\Lambda>0$ is zero as we expected. Furthermore, for
$\Lambda=0$, this solution gives the pressure of a universe filled by a
perfect fluid. We call the solutions (5.7), (5.8), (5.9) and (5.13) in this
article as cosmological solutions of the BD equations. It should be noted that
these solutions were previously obtained in less detail without the use of the
concept of symmetry by one of the authors of this article in [41, 42]. Note
that when the coupling constant $\omega$ tends to infinity in solutions (5.7),
(5.8), (5.9) and (5.13), one concludes that they are actually the cosmological
solutions of general relativity, which can be obtained from solving EFEs in
the presence of the cosmological constant for a universe filled with a perfect
fluid with the equation of state $p=w\rho$.

Before closing this section, let us assume that the perfect fluid, filling the
universe, is radiation. For $w=1/3$ and $\Lambda=0$ it follows from (5.6) that
$n=0$. Then, equations (3.19) and (3.20) are reduced to

$\displaystyle\frac{8\pi\rho}{C}$ $\displaystyle=$ $\displaystyle
3\frac{\dot{a}^{2}}{a^{2}},$ (5.14) $\displaystyle\frac{8\pi p}{C}$
$\displaystyle=$
$\displaystyle-2\frac{\ddot{a}}{a}-\frac{\dot{a}^{2}}{a^{2}}.$ (5.15)

Now we compare these equations with the case $\Lambda=0$ of the following
Friedmann equations in the theory of general relativity

$\displaystyle 8\pi G\rho$ $\displaystyle=$ $\displaystyle
3\frac{\dot{a}^{2}}{a^{2}},$ (5.16) $\displaystyle 8\pi G\rho$
$\displaystyle=$
$\displaystyle-2\frac{\ddot{a}}{a}-\frac{\dot{a}^{2}}{a^{2}}.$ (5.17)

In order to be compatible with general relativity in this special case (for
radiation in the absence of the $\Lambda$), the value of the constant $C$
should be

$\displaystyle C(w=\frac{1}{3},\omega,G,\Lambda=0)=\frac{1}{G}.$ (5.18)

By using condition (5.18), one can obtain the cosmological solutions of the BD
equations for radiation in the absence of the $\Lambda$, giving us

$\displaystyle
a(t)={c^{\prime}}_{1}\sqrt{2(t+{c^{\prime}}_{2})},~{}~{}~{}~{}\varphi(t)=\frac{1}{G},~{}~{}~{}~{}p(t)=\frac{1}{32\pi
Gt^{2}},~{}~{}~{}~{}~{}~{}~{}~{}\rho(t)={3}p(t).$ (5.19)

Therefore, when the condition (5.18) holds, then not only general relativity
and Brans-Dicke theories are compatible with each other, but also the
solutions of both theories will be the same for radiation. This result is
indicative of the fact that for a universe full of radiation $w=1/3$ whose
metric is FLRW flat $(k=0)$, in the absence of the cosmological constant
$\Lambda$, the coupling constant of the BD scalar field $\varphi$ with the
gravitational field $g_{\mu\nu}$ is at its lowest possible, so that the two
fields $\varphi$ and $g_{\mu\nu}$ are completely separate and have no paining
with other. In other words, the curvature of space-time is very large compared
to the curvature of the scalar field.

##  6 Special solutions of the FBD-equation

In order to make sure the correctness of cosmological solutions of the BD
equations (5.7), (5.8), (5.9) and (5.13), it is necessary to compare our
solutions with special solutions such as Nariai’s solutions [43], O’Hanlon-
Tupper vacuum solutions [8] and the inflation solutions that have been
obtained before. It can be shown that when $w\in[0,1/3]$ and
${c^{\prime}}_{2}=0$, our solutions include the Nariai’s solutions. Nariai’s
solutions in the absence of the $\Lambda$, when the fluid constituting the
universe is dust, include the BD dust solutions

$\displaystyle a(t)$ $\displaystyle=$
$\displaystyle{c^{\prime}}_{1}\Big{(}\frac{4+3\omega}{2(\omega+1)}\Big{)}^{\frac{2(\omega+1)}{4+3\omega}}~{}~{}(t+{c^{\prime}}_{2})^{\frac{2(\omega+1)}{4+3\omega}},$
(6.1) $\displaystyle\varphi(t)$ $\displaystyle=$ $\displaystyle
C\Big{(}{c^{\prime}}_{1}\frac{(4+3\omega)}{2(\omega+1)}\Big{)}^{\frac{2}{4+3\omega}}~{}~{}(t+{c^{\prime}}_{2})^{\frac{2}{4+3\omega}},$
(6.2) $\displaystyle\rho(t)$ $\displaystyle=$ $\displaystyle
M(\omega,0)~{}X(\omega,0)~{}~{}~{}(t+{c^{\prime}}_{2})^{\frac{-6(\omega+1)}{4+3\omega}},$
(6.3) $\displaystyle p(t)$ $\displaystyle=$ $\displaystyle 0.$ (6.4)

where $\omega\neq-4/3,-1.$ To obtain the vacuum solution ($p=0$, $\rho=0$), we
reconsider equation (3.19) that for the vanishing energy density we have

$\displaystyle 3-3n-\frac{\omega}{2}n^{2}=0,$ (6.5)

The solutions of this quadratic equation are:

$\displaystyle n_{\pm}=\frac{-3\pm\sqrt{3(3+2\omega)}}{\omega}.$ (6.6)

The solutions $a(t)$ and $\phi(t)$ obtained for these two values of $n$ are
called the “vacuum solutions” of the BD equations. By substituting (5.6) into
equation (6.6) one may obtain the parameter of the equation of state $w$ in
terms of the coupling parameter $\omega$, expressing

$\displaystyle
w_{\pm}=\frac{-(2\omega+3)\pm(\omega+1)\sqrt{3(3+2\omega)}}{\pm\omega\sqrt{3(3+2\omega)}}.$
(6.7)

For these two values of $w_{\pm}$, the cosmological solutions of BD equations
in the absence of the $\Lambda$ are worked out

$\displaystyle a_{\pm}(t)$ $\displaystyle=$ $\displaystyle
a_{0}(1+bt)^{q_{\pm}},$ (6.8) $\displaystyle\varphi_{\pm}(t)$ $\displaystyle=$
$\displaystyle\varphi_{0}(1+bt)^{1-3q_{\pm}},$ (6.9)

where

$\displaystyle
q_{\pm}=\frac{\omega+1\pm\sqrt{\frac{2\omega+3}{3}}}{3\omega+4},$

also, the constants $a_{0}$, $\varphi_{0}$ and $b$ are

$\displaystyle
a_{0}:=a(t=0)={c^{\prime}}_{1}\Big{(}\frac{{c^{\prime}}_{2}}{q_{\pm}}\Big{)}^{q_{\pm}},~{}~{}~{}~{}~{}~{}\varphi_{0}:=\varphi(t=0)=\frac{1}{G}~{}{{c^{\prime}}_{1}}^{\frac{1-3q_{\pm}}{q_{\pm}}}\Big{(}\frac{{c^{\prime}}_{2}}{q_{\pm}}\Big{)}^{1-3q_{\pm}},~{}~{}~{}~{}~{}~{}b:=\frac{1}{{c^{\prime}}_{2}},$

we note that the range of parameter $\omega$ is $\omega>-3/2$ with
$\omega\neq-4/3,0$. These solutions are in agreement with the solutions of
those of [15]. The special solutions (6.8) and (6.9) were first obtained by
O’Hanlon and Tupper in 1972 and are usually called the vacuum solutions of
O’Hanlon-Tupper [8]. These solutions are completely consistent with Chauvet’s
solutions [9].

At the end of this section, let us examine the inflation solutions of BD
equations. For the state parameter $w=-1$, the universe was passing through
the false vacuum era [44, 45, 46]. We expect that the solutions obtained from
the BD equations in this particular state to be inflationary. To this end, we
begin with the case $\Lambda=0$ of the cosmological solutions (5.7), (5.8) and
(5.9) for $w=-1$, which are, respectively, given by

$\displaystyle a(t)$ $\displaystyle=$
$\displaystyle{c^{\prime}}_{1}\Big{(}\frac{1}{\omega+\frac{1}{2}}\Big{)}^{\omega+\frac{1}{2}}~{}~{}(t+{c^{\prime}}_{2})^{\omega+\frac{1}{2}},$
(6.10) $\displaystyle\varphi(t)$ $\displaystyle=$
$\displaystyle\frac{1}{G}~{}{c^{\prime}}_{1}^{\frac{2}{\omega+\frac{1}{2}}}\Big{(}\frac{2}{2\omega+1}\Big{)}^{2}~{}~{}(t+{c^{\prime}}_{2})^{2},$
(6.11) $\displaystyle\rho(t)$ $\displaystyle=$ $\displaystyle\frac{1}{8\pi
G}~{}{c^{\prime}}_{1}^{\frac{2}{\omega+\frac{1}{2}}}~{}~{}\frac{(2\omega+3)(6\omega+5)}{(2\omega+1)^{2}}:=\rho_{-1}.$
(6.12)

where $\omega\neq-1/2,-3/2,-5/6$. The last one shows that the energy density
of the universe in the false vacuum era is a constant. Here, we have denoted
this constant value by $\rho_{-1}$. Now, it is necessary to express the
integration constant ${c^{\prime}}_{1}$ in equations (6.10) and (6.11) in
terms of cosmic scale factor of the universe at the time of the Big Bang
$t=0$, in which we denote by $a_{0}$. For this purpose, by putting $t=0$ in
cosmological solution (6.10) and according to the initial condition
$a_{0}=a(t=0)$, we then get

$\displaystyle{c^{\prime}}_{2}=\frac{2\omega+1}{2}\Big{(}\frac{a_{0}}{{c^{\prime}}_{1}}\Big{)}^{\frac{2}{2\omega+1}}.$
(6.13)

By substituting (6.12) and (6.13) into solutions (6.10) and (6.11), one can
obtain

$\displaystyle a(t)$ $\displaystyle=$ $\displaystyle a_{0}(1+\chi
t)^{\omega+\frac{1}{2}},$ (6.14) $\displaystyle\varphi(t)$ $\displaystyle=$
$\displaystyle\frac{1}{G}~{}a_{0}^{\frac{2}{\omega+\frac{1}{2}}}~{}(1+\chi
t)^{2},$ (6.15)

where

$\displaystyle\chi^{2}:=\frac{32\pi
G\rho_{-1}~{}{c^{\prime}}_{1}^{-\frac{2}{\omega+\frac{1}{2}}}}{(6\omega+5)(2\omega+3)},$
(6.16)

is a constant. The cosmic scale factor (6.14) and the BD scalar field (6.15)
are the inflation solutions of the BD equations in the absence of the
cosmological constant, respectively. When the coupling parameter $\omega$
tends to infinity, for very small times $\chi t\ll 1$, the limits of the
inflation solutions of the BD equations become

$\displaystyle\lim_{{}_{\omega\rightarrow\infty}}a(t)$ $\displaystyle=$
$\displaystyle
a_{0}~{}\lim_{{}_{\omega\rightarrow\infty}}\Big{(}1+\frac{1}{\omega}\sqrt{\frac{8\pi
G\rho_{-1}}{3}}t\Big{)}^{\omega}=a_{0}~{}e^{\sqrt{\frac{8\pi
G\rho_{-1}}{3}}t},$ (6.17)
$\displaystyle\lim_{{}_{\omega\rightarrow\infty}}\varphi(t)$ $\displaystyle=$
$\displaystyle\frac{1}{G}.$ (6.18)

It can be seen that in the limit $\omega\rightarrow\infty$, the cosmic scale
factor $a(t)$ changes exponentially with cosmic time $t$, while the BD scalar
field $\varphi(t)$ is a constant. Notice that solutions (6.17) and (6.18) are
the same solutions which are obtained from the theory of the GR in the absence
of the cosmological constant in the case where the equation of state of the
universe is in the form $p=-\rho_{-1}$. For the case that $\chi t\gg 1$, then
the inflationary solutions (6.14) and (6.15) are as follows:

$\displaystyle a(t)$ $\displaystyle=$ $\displaystyle a_{0}~{}(\chi
t)^{\omega+\frac{1}{2}},$ (6.19) $\displaystyle\varphi(t)$ $\displaystyle=$
$\displaystyle\varphi_{0}~{}(\chi t)^{2},$ (6.20)

where

$\displaystyle\varphi_{0}=\frac{1}{G}~{}a_{0}^{\frac{2}{\omega+\frac{1}{2}}}.$
(6.21)

For the universe in the false vacuum era $(p=-\rho_{{}_{-1}})$ and in the
absence of cosmological constant $\Lambda$, it follows from (5.6) that

$\displaystyle n_{|_{w=-1,\Lambda=0}}=-\frac{2}{\omega+\frac{1}{2}}.$ (6.22)

Then, by substituting (6.22) into the power-law equation (3.18) and by using
the fact that $C=1/G$, we find

$\displaystyle\frac{1}{G}=C=\varphi(t)~{}a^{n}(t)=\varphi(t=0)~{}a^{n}(t=0)=\varphi_{0}~{}a_{0}^{n}=\varphi_{0}a_{0}^{-\frac{2}{\omega+\frac{1}{2}}}.$

which is nothing but equation (6.21). This issue can be a confirmation of the
correctness of the inflationary solutions obtained by solving the BD equations
for the false vacuum universe in the absence of the cosmological constant.

##  7 Conclusions

In this work, we have focused on solving the BD equations for spatially flat
$(k=0)$ FLRW universe fulfilling a perfect fluid with the equation of state
$p=\omega\rho$, $(-1\leq\omega\leq 1)$. Using equations (3.12) and (3.18) we
have solved the BD equations in both cases of the absence and presence of the
cosmological constant to obtain functions $a(t)$, $\varphi(t)$, $\rho(t)$ and
$p(t)$. As was shown, the constant $n$ in equation (3.18) and the state
equation parameter $w$ in (3.12) are compatible with each other, provided that
the condition (5.6) is held. Equations (3.12) and (3.18) helped us summarize
all BD equations in one equation, (3.25), which we called FBD-equation. This
equation alone could describe the dynamics of the spatially flat FLRW
cosmological model in the framework of BD theory. As we have seen, this
cosmological model could be imagined as a $1$-dimensional dynamical system
with the configuration space $\mathbb{Q}=(a)$ in which the governing equation
of particle motion was given by
$\mathbf{F}=\phi(t,a,\dot{a})\mathbf{\partial}_{a}$. As the first step to
solve the FBD-equation we found the corresponding Lie point symmetries. Then,
we showed that the FBD-equation has two independent point symmetries with
infinitesimal generators $\mathbf{X}_{1}=\partial_{t}$ and
$\mathbf{X}_{2}=a\partial_{a}$. Using DPs and $\lambda$-symmetries, we found
two independent $2$-tuple $(S_{1},R_{1})$ (equation (4.7)) and $(S_{2},R_{2})$
(equation (4.9)) consisting of null form $S$ and integrating factor $R$. Using
these independent $2$-tuples and applying Duarte’s integral formula (2.4), we
obtained the two independent invariants (4.8) and (4.10) which were associated
to the infinitesimal generators $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$,
respectively. We showed that these two independent invariants, each of which
was a first-order ODE for the cosmic scale factor $a(t)$ in terms of cosmic
time $t$, could be viewed as a system of algebraic equations for the unknown
functions $a(t)$ and $\dot{a}(t)$. By eliminating $\dot{a}(t)$ between the two
equations of the system, we obtained the analytical solution (5.7) for the
cosmological scale factor $a(t)$. By employing this solution in the BD
equations, we were able to obtain the BD scalar field $\varphi(t)$, the energy
density of the universe $\rho(t)$ and its pressure $p(t)$, which were given by
equations (5.8), (5.9), and (5.13), respectively. Not only these solutions
gave the scale factor $a(t)$, the BD scalar field $\varphi(t)$, the energy
density $\rho(t)$ and pressure of the universe $p(t)$ for the dust, but also
they gave us the more general state of a perfect fluid with equation of state
$p=\omega\rho$, $(-1\leq\omega\leq 1)$.

As an interesting result, we have shown that when the coupling parameter
$\omega$ tends to infinity, the cosmological solutions of GR theory can be
obtained from the analytical solutions of BD equations (5.7), (5.8), (5.9) and
(5.13); this can be a confirmation of $\lim_{\omega\rightarrow\infty}$ BD =
GR. Note that our cosmological solutions are rich, so that they include many
well-known special solutions that were previously obtained by other methods.
For example, the solutions presented by Uehara and Kim [12] are special cases
of our solutions, when the cosmological constant $\Lambda$ is present for the
dust dominated universe. In addition, as shown, the case $\Lambda=0$ of our
solutions include the following famous solutions: Nariai’s solutions
$(w\in[0,1/3])$, BD dust solutions $(w=0)$, radiation solutions $(w=1/3)$,
inflationary solutions $(w=-1)$, O’Hanlon-Tupper vacuum solutions
$(w=w_{\pm})$ and general relativity solutions $(\omega\rightarrow\infty)$.

### Acknowledgements

This work has been supported by the research vice chancellor of Azarbaijan
Shahid Madani University under research fund No. 1402/917.

## References

  * [1] J. B. Griffiths, J. Podolsky, Exact space-times in Einstein’s general relativity, New York: Cambridge University Press, (2009); E. P. Wigner, Events, laws of nature, and invariants principles, Science, 145 (1964) 995-999.
  * [2] Y. Fujii, and K. Maeda, The scalar-tensor theory of gravitation, Cambridge: Cambridge University Press, (2007).
  * [3] C. Brans and R. H. Dicke, Mach’s principle and a relativistic theory of gravitation, Phys. Rev. (series I), 124(3) (1961) 925-935.
  * [4] C. H. Brans, Mach’s principle and a relativistic theory of gravitation. II, Phys. Rev.  (series I), 125(6) (1962) 2194-2201.
  * [5] P. A. M. Dirac, Cosmological models and the large numbers hypothesis, Proceedings of the Royal Society of London A, 338 (1974) 439-446.
  * [6] A. B. Batista, Inhomogeneous distributions in the Brans-Dicke cosmology, Phys. Rev. D, 21(8) (1980) 2119-2121.
  * [7] C. B. G. MacIntosh, Position dependent Robertson-Walker solutions of the Brans-Dicke field equations, Phys. Lett. A, 43(1) (1973) 33-34.
  * [8] J. O’Hanlon, B. O. J. Tupper, Vacuum-field solutions in the Brans-Dicke theory, IL Nuovo Cimento B, (1971-1996), 7(2) (1972) 305-312.
  * [9] P. Chauvet, General vacuum-universe solutions in Brans-Dicke cosmology, Astrophys. Space Sci., 90(1) (1983) 51-58.
  * [10] J. M. Cervero and P. G. Estevez, General solutions for a cosmological Robertson-Walker metric in the Brans-Dicke theory, Gen. Relativ. Gravit.,  15(4) (1983) 351-356.
  * [11] D. Lorentz-Petzold, Comment on the general vacuum solutions in Brans-Dicke cosmology, Astrophys.Space Sci., 96(2) (1983) 451-453.
  * [12] K. Uehara and C. W. Kim, Brans-Dicke cosmology with cosmological constant, Phys. Rev. D, 26(10) (1982) 2575-2579.
  * [13] S. M. M. Rasouli and P. V. Moniz, Exact cosmological solutions in modified Brans-Dicke theory, Class. Quantum Gravity, 33(3) (2016) 035006.
  * [14] K. Atazadeh and F. Darabi, $f(T)$ cosmology via Noether symmetry, Eur. Phys. J. C, 72(5) (2012) 1-7.
  * [15] S. Capozziello and V. Faraoni, Beyond Einstein Gravity: A survey of gravitational theories for cosmology and astrophysics, Heidelberg: Springer Science+Business Media, (2011).
  * [16] P. Rastall, Generalization of the Einstein theory, Phys. Rev. D, 6 (1972) 3357.
  * [17] V. K. Chandrasekar, M. Senthilvelan and M. Lakshmanan, On the complete integrability and linearization of certain second-order nonlinear ordinary differential equations, Proc. R. Soc. A, 461 (2005) 2451-2476.
  * [18] V. K. Chandrasekar, M. Senthilvelan and M. Lakshmanan, Extended Prelle-Singer method and integrability/solvability of a class of nonlinear nth order ordinary differential equations, J. Non lin. Math. Phys., 12 (suppl.) (2005) 184-201.
  * [19] R. Mohanasubha, V. K. Chandrasekar, M. Senthilvelan and M. Lakshmanan, Interplay of symmetries, null forms null, Darboux polynomials, integrating factors and Jacobi multipliers in integrable second-order differential equations, Proc. R. Soc. A, 470(2163) (2013) 20130656.
  * [20] R. Mohanasubha, V. K. Chandrasekar, M. Senthilvelan, and M. Lakshmanan, Interconnections between various analytic analytic approaches applicable to third-order nonlinear differential equations, Proc. R. Soc. A, 471(2176) (2015) 20140720.
  * [21] M. Manoranjani, R. Mohanasubha, V. K. Chandrasekar and M. Senthilvelan, Extended Prelle-Singer procedure and Darboux polynomials method: An unknown interconnection, Int. J. Non Lin. Mech., 118 (2020) 103284.
  * [22] V. K. Chandrasekar, M. Senthilvelan and M. Lakshmanan, M., On the complete integrability and lineariz- ation of certain second-order nonlinear ordinary differential equations, Proc. R. Soc. A, 461(2060) (2005) 2451-2476.
  * [23] C. Muriel, and J. L. Romero, New methods of reduction for ordinary differential equations, IMA J. Appl. Math., 66(2) (2001) 111-125.
  * [24] C. Muriel and J. L. Romero, First integrals, integrating factors, and $\lambda$-symmetries of second-order Differ- ential equations, J. Phys. A: Math. Theor., 42(36) (2009) 365207.
  * [25] C. Muriel and J. L. Romero, Nonlocal symmetries, telescopic vector fields and $\lambda$-symmetries of ordinary differential equations, SIGMA, 8 (2012) 106-127.
  * [26] A. Bhuvaneswari, R. A. Kraenkel and M. Senthilvelan, Application of the $\lambda$-symmetries approach and time independent integral of the modified Emden equation, Nonlin. Anal.: Real World Appl, 13(3) (2012) 1102-1114.
  * [27] E. Ahmadi-Azar, K. Atazadeh and A. Eghbali, Exact cosmological solutions via invariants of the symmetry groups, arXiv:2310.19772 [gr-qc].
  * [28] L. G. S. Duarte, S. E. S. Duarte, A. C. P. Da Mota and J. E. F. Skea, Solving the second-order ordinary differential equations by extending the Prelle-Singer method, J. Phys. A: Math. Gen.,  34(14) (2001) 3015-3024.
  * [29] G. W. Bluman and S. Komei, Symmetries and differential equations, Vol. 81, New York: Springer Science+Business Media, (1989).
  * [30] N. H. Ibragimov, Elementary Lie groups and ordinary differential equations, Chichester: John Wiley, (1999).
  * [31] S. Capozziello, M. De Laurentis, Extended theories of gravity, Phys. Rep., 509(4-5) (2011) 167-321.
  * [32] C. M. Will, The confrontation between general relativity and experiment: An update, Phys. Rep., 113(6) (1984) 345-422.
  * [33] V. Faraoni, The $\omega\rightarrow\infty$ limit of Brans-Dicke theory, Phys. Lett. A, 245 (1998) 26-30.
  * [34] N. Banerjee and S. Sen, Does Brans-Dicke theory always yield general relativity in the infinite $\omega$ limit? Phys. Rev. D, 56(2) (1997) 1334-1337.
  * [35] H. Dehnen, and O. Obrego’s, Ecact cosmological solutions in Brans-Dicke’s scalar-tensor theory, I, Astrophys. Space Sci., 14(2) (1971) 454-459.
  * [36] E. Gourgoulhon, 3+1 formalism in general relativity, Springer, (2012).
  * [37] P. Chvanis, Cosmology with a stiff matter era, Phys. Rev. D, 92(10) (2015) 103004-103034.
  * [38] G. Brando, J. C. Fabris, F. T. Falciano and O. Galkina, Stiff matter solution in Brans-Dicke theory and the general relativity limit, Int. J. Mod. Phys. D, 28(12) (2019) 1950156.
  * [39] V. B. Johri and D. Kalyyani, Cosmological models with constant deceleration parameter in Brans-Dicke theory, Gen. Rela. Grav., 26(12) (1994) 1217-1232.
  * [40] H. Kim, Brans-Dicke as a unified model for dark matter-dark energy, Mon. Not. R. Astro. Soc., 364(3) (2005) 813-822.
  * [41] E. Ahmadi-Azar and N. Riazi, A class of cosmological solutions of Brans-Dicke theory with cosmological constant, Astrophys. Space Sci., 226(1) (1995) 1-5.
  * [42] I. Ahmadi-Azar, Cosmological solutions of Brans-Dicke equations with cosmological constant, Iranian J. Phys. Research, 3(3) (2002) 177-183.
  * [43] H. Nariai, On the Green’s function in an expanding universe and its role in the problem of Mach’s principle, Prog. Theo. Phys., 40(1) (1968) 49-59.
  * [44] C. Mathiazhagan and V. B. Johri, An inflationary universe in Brans-Dicke theory: a hopeful sign of theortetical estimation of the gravitational constant, Class. Quantum Gravity, 1 (1984) L29.
  * [45] J. D. Barrow and K. Maeda, Extended inflationary universes, Nucl. Phys. B, 341 (1990) 294-308.
  * [46] R. G. Vishwakarma, A Machian model of dark energy, Class. Quantum Gravity,  19(18) (2002) 4747-4752.


MLT-LE: predicting drug–target binding affinity with multi-task residual
neural networks

Elizaveta Vinogradova1,2, Karina Pats1,2, Ferdinand Molnár2*, Siamac Fazli3*,

1 Department of Biology, Nazarbayev University, Nur-Sultan, Kazakhstan
2 Computer Technologies Laboratory, ITMO University, Russia
3 Department of Computer Science, Nazarbayev University, Nur-Sultan,
Kazakhstan

*<EMAIL_ADDRESS><EMAIL_ADDRESS>
## Abstract

Assessing drug-target affinity is a critical step in the drug discovery and
development process, but to obtain such data experimentally is both time
consuming and expensive. For this reason, computational methods for predicting
binding strength are being widely developed. However, these methods typically
use a single-task approach for prediction, thus ignoring the additional
information that can be extracted from the data and used to drive the learning
process. Thereafter in this work, we present a multi-task approach for binding
strength prediction. Our results suggest that these prediction can indeed
benefit from a multi-task learning approach, by utilizing added information
from related tasks and multi-task induced regularization.

Availability of implementation: Associated data, pre-trained models, and
source code are publicly available at https://github.com/VeaLi/MLT-LE.

## Introduction

One of the early steps in drug discovery is the identification of compounds
that have a strong binding to a given target protein. However, due to high
costs and time constraints, it is not possible to screen a large subset of the
chemical space by performing laboratory experiments. To solve this problem,
many computational methods based on machine learning have been developed [1,
2, 3, 4, 5, 6, 7, 8, 9]. Most of them are designed to predict one binding
constant at a time [1, 3, 5, 6]. However, these methods may struggle to adapt
to a new drug-target affinity (DTA) prediction pair due to inconsistency and
under-representation of training data [2, 10, 11]. To alleviate this problem,
multi-task learning methods have recently been employed [2].

The multi-task approach to molecular property prediction takes advantage of
all consistent data by predicting multiple properties simultaneously. By
utilizing multiple binding constants at once, multi-task prediction models
show improved ability to adapt to new data [2, 10, 12, 13, 14, 7, 15]. More
technically, multi-task learning is an approach in machine learning that
forces a model to favor certain hypotheses over others [12, 16]. Specifically,
multi-task learning biases the model toward hypotheses that contribute to the
solution of more than one task. Thus, the model’s attention becomes focused on
the features that really matter, since multiple tasks provide additional
evidence for their relevance or irrelevance, thus ideally providing additional
regularization. Therefore, in general, models trained to solve multiple tasks
are more likely to escape the data-dependent noise and find more useful
representations [16]. In addition, models trained to solve multiple related
tasks simultaneously show less tendency to overfit, generalize better, and
have access to more sample data [13, 11, 14, 16, 17].

Multi-task learning can be seen as being inspired by human learning. Humans
are able to learn sequentially without forgetting previous knowledge and tend
to use the knowledge they already acquired to solve different new tasks [16,
18]. Initially, multi-task learning was introduced to regularize parameters in
linear models and evaluate the task relatedness [16]. Nowadays, multi-task
learning is used to find solutions to complex problems involving combinations
of multiple data sources, such as in bioinformatics [19]. For example, in
[19], it was shown that multi-task learning can improve the quality of gene
regulatory network reconstruction in settings of scarcely labeled data by
using data from multiple organisms simultaneously relying on orthologous
genes.

The protein-ligand binding strength can be measured in terms of different
binding constants. As a result, there exist several independent datasets.
Generally, these datasets are used to train different single-task algorithms.
Combining this data would naturally produce more powerful prediction models
[12]. However, only a few multi-task drug target binding strength prediction
methods have been developed to date [2, 7]. Why is this? Mainly because only
few drug target pairs have several measurements available in terms of binding
constants. For example, for a dataset with $\sim 717$ thousand records taken
from BindingDB, the complete overlap between the four most common binding
constants $\\{K_{d},K_{i},IC_{50},EC_{50}\\}$ is only $96$ records (Fig. 1).

Figure 1:  Figure 1. The overlap across all four most commonly expressed
binding constants for human drug-target interaction data found in BindingDB is
only 96 records (center).

Despite this, a few multi-task methods for drug-protein binding prediction [2,
7] have shown that they are able to exploit multiple overlaps (intersection)
of the data and outperform single-task models in multiple ways [2, 7].

In this work, we propose a new multi-task framework to address the challenges
encountered in these related studies [2, 7]. In particular, we propose an
approach that allows using all available data (union), including missing
values, and using auxiliary tasks to improve performance. We compare our
approach with the state-of-the-art GraphDTA framework, which has been shown to
outperform most existing methods [5].

## Materials and Methods

### Overview of MLT-LE

We propose a new deep learning framework - MLT-LE for the prediction of
binding strength. Existing methods render the binding strength prediction task
as a single-task problem, whereas we propose to solve it in a multi-task
manner (Fig. 2). The models in our framework are highly adjustable and support
multiple drug and target encodings.

Figure 2:  General structure of MLT-LE model. The basic structure of the model
within MLT-LE is implemented as a residual convolutional network with
adjustable depth and number of outputs.

### Handling of missing data

Combining human data from the BindingDB database [20] for any subset of
binding constants almost always results in $>50$% missing data (see table 3
for an example), which makes multi-task approaches challenging to implement.
Especially since it has to be solved during algorithm training rather than
during the preparation phase, since removing the missing data would greatly
reduce the sample size. At the same time, filling in the missing data on such
a large scale is impossible. Dahl _et al._ suggested to simply not propagate
the error for missing values [17] and Hu _et al._ proposed to use branching -
propagating the error only for non-missing values using a separate branch in
the neural network [2]. However, usually, when encountering the problem of
training an algorithm with missing data points, a technique called masking is
used [17]. Masking is the process of dynamically hiding the missing data (i.e.
excluding it from the computations) whenever the error is calculated (which is
well described for batch optimizers). The error is calculated using some loss
function, so masking is part of the loss calculation. In this work, a
different loss adjustment is suggested that shows promising results compared
to the masking technique. The developed approach proposes to propagate the
last non-missing discounted error found, if there is no data available, thus
allowing the use of stochastic optimizers.

### Benchmark

To assess the performance of our models, we compare them to the state-of-the-
art GraphDTA framework [5]. We used the model hyperparameters for the GraphDTA
models (Table 1), as suggested in the original work [5], and compared the
performance of all models on the same split of training, validation and test
data. Training was conducted for $1000$ epochs and the best models were
selected based on the performance on the validation data, particularly by
monitoring changes in the mean squared error (MSE) for the target variable(s).
Subsequently, all the best performing models were used to obtain predictions
for unseen sets of test data.

Table 1: Models hyperparameters. Framework |  Optimizer | Learning rate | Batch size
---|---|---|---
GraphDTA |  Adam | 0.0005 | 1024
MLT-LE |  LookAhead/Nadam with sync. period of 3 | 0.001 | 1024

### Complexity

Table 2 gives an overview of the complexity of the models involved in the
experiments. As can be seen, the total number of parameters for all models is
comparable, although some models solve more tasks. The numbers in the model
names indicate the number of layers. The plus sign in the model names
indicates that the model is trained on merged $K_{d}$+$EC_{50}$ data (see
Section ”Data”). The prefix ’Res’ indicates that the network has residual
connections. For more details on the GraphDTA architecture, please refer to
[5].

Table 2: Model description. Framework |  Model name |  Drug representation |  Protein representation |  Num. of trainable parameters |  Num. of tasks
---|---|---|---|---|---
GraphDTA |  GAT2 |  Graph attention |  Convolution |  1,460,845 |  1
GATGCN2 |  Graph attention - graph convolution combined |  |  4,749,573
GCN3 |  Graph convolution |  |  2,062,767
GIN5 |  Graph isomorphism |  |  1,297,505
MLT-LE |  ResCNN1 |  Convolution |  Convolution |  2,391,783 |  7
ResCNN1+ |  |  |
ResCNN1GCN4 |  Graph convolution |  |  2,627,895
ResCNN1GCN4+ |  |  |
ResCNN1GIN5 |  Graph isomorphism |  |  2,290,857
ResCNN1GIN5+ |  |  |

### Data

In this work, we used human $K_{d}$ and $EC_{50}$ data from BindingDB v2022m3
[20], each divided into $3$ sets of the same size, namely train, validation,
and test, to mitigate the effects of overfitting. The original subset of human
data from BindingDB contained $1,274,797$ records, with $525,152$ unique
SMILES and $1,981$ unique amino acid sequences of target proteins. Only
records with at least one binding constant were kept. Then, the records were
filtered to contain only binding constants recorded without range boundaries,
indicated by $+-$ and $><$ signs. The SMILES sequences were then canonicalized
using RDkit-pypi v2022.3.5 [21], and the isomeric information was removed.
Only records with valid SMILES strings were kept. As a result, the dataset
contained slightly fewer unique SMILES sequences, $498,742$. The QED
coefficient [22] was then calculated for each SMILES using RDkit (to be used
as an auxiliary task later). The individual data records were then aggregated
using the median, as was done in [10]. Only valid protein sequences starting
with a start codon were considered. Then, the individual records were labeled
as “active” or “non-active” drug-target pairs using a threshold of
$(K_{d}|K_{i}|IC_{50}|EC_{50})<1\mu M$ for binarization, as done in [23, 24]
or similarly in [25]. After that, the training dataset contained a total of
$716,342$ records. The overlap of records containing different binding
constants was then calculated (Fig. 1); the total overlap between records was
$96$ records in total, with missing values corresponding to more than $>77\%$.
Next, only records with $K_{d}$ or $EC_{50}$ non-missing values were kept,
yielding $19,796$ $K_{d}$ records and $63,487$ $EC_{50}$ records. The data
were then shuffled and divided into three equal parts for each set, resulting
in $6,598$ training, $6,598$ validation, and $6,600$ test records for $K_{d}$
and $21,153$, $21,153$, $21,155$ records for $EC_{50}$ respectively, all of
which contain non-overlapping drug-target pairs. In addition, we created a
combined $K_{d}$+$EC_{50}$ dataset, denoted in this work by ”+”, by combining
the training and validation datasets from the separate $K_{d}$ and $EC_{50}$
data. It should be noted that we also saved the target pH information for all
datasets in order to use it as another auxiliary task. Table 3 shows the
number of non-missing values in the train data for each dataset. Note that the
validation and test datasets have a similar distribution of records as the
corresponding train dataset.

Table 3: Number of non-missing values in train data. Dataset |  $\mathbf{K_{d}}$ |  $\mathbf{K_{i}}$ |  $\mathbf{IC_{50}}$ |  $\mathbf{EC_{50}}$ |  active |  pH |  QED
---|---|---|---|---|---|---|---
$K_{d}$ set |  6,598 |  11 |  20 |  1 |  6,598 |  580 |  6,598
$EC_{50}$ set |  3 |  159 |  201 |  21,153 |  21,153 |  1,144 |  21,153
Combined (’$+$’) $K_{d}+EC_{50}$ set |  6,600 |  170 |  221 |  21,153 |  27,750 |  1,723 |  27,750

## Results

Figure 3 shows the performance of the best models obtained on the test set. It
can be seen that the MLT-LE models showed comparable or better performance
than the baseline GraphDTA models on all metrics. Additional statistical
testing of the observed difference between the predictions of the different
frameworks can be found in the supplementary materials. The best scores
obtained for each metric are shown in bold.

Figure 3:  Comparison of the model performances on the test set for $K_{d}$
and $EC_{50}$ data. The best scores obtained for each metric are shown in
bold.

Figure 4 shows the performance of all models used in the experiment on the
validation data. It can be seen that the GraphDTA and MLT-LE framework models
show comparable performance, although the GraphDTA models show more stable
convergence, while MLT-LE achieves both the highest concordance index ($CI$)
and solves more tasks.

It was also observed that the performance of MLT-LE models is also of high
quality for other tasks; more information on this is provided in the
supplementary material.

In addition, similar model performance was observed when training with a
smaller batch size of $256$ over $100$ epochs for the same datasets (the
results of this experiment are shown in the supplementary material).

Figure 4:  Comparison of model performance on the validation data for $K_{d}$
and $EC_{50}$ data. GraphDTA performance is shown in orange and MLT-LE
performance is shown in blue.

## Discussion

Incorporating data from various sources helps to obtain more accurate
predictions of molecular properties [12]. However for this task, the
combination of data sources results in a substantial amount of missing data
(as can be seen in Table 3). Thus, the question may arise whether combining
these data would be beneficial for the prediction [7, 15]. The results of the
present work suggest that missing data can be included and is shown to be
beneficial in terms of performance (Fig. 3). We attribute this observation to
the effect of multi-task models being able to incorporate data from multiple
sources and thus utilizing information that is not available for single-task
models.

## Conclusion

In this work, we compared our binding strength prediction framework - MLT-LE
with the state-of-the-art GraphDTA framework. We tested the performance on two
different subsets of BindingDB. We found that MLT-LE performs well for the
benchmark datasets across all metrics. We furthermore observed that a multi-
task convolutional model can provide comparable results to a single-task graph
model, and that a multi-task graph model can perform better than a single-task
graph model. We also observed the beneficial effects of incorporating more
data, even given the high percentage of missing values.

We believe that it may be possible to further improve performance on the
binding strength prediction task using data from more organisms if the
orthologous protein sequences in these organisms are matched properly and if a
special task is used in the neural network to distinguish the organisms, which
can smooth out the effects of different medians of binding constants in these
organisms.

## Supporting Information

Associated and supplementary data, data preparation routines, pre-trained
models, and source code are publicly available at
https://github.com/VeaLi/MLT-LE.

## Acknowledgments

The co-authors would like to acknowledge the support of ITMO University
Research Grants funding to EV and KP (621314) and Nazarbayev University
Research Grants funding to SF and FM (240919FD3926, 091019CRP2108 and
110119FD4520).

## References

  * 1. T. He, M. Heidemeyer, F. Ban, A. Cherkasov, and M. Ester, “SimBoost: a read-across approach for predicting drug–target binding affinities using gradient boosting machines,” _J. Cheminformatics_ , vol. 9, no. 1, 2017.
  * 2. F. Hu, J. Jiang, D. Wang, M. Zhu, and P. Yin, “Multi-PLI: interpretable multi‐task deep learning model for unifying protein–ligand interaction datasets,” _J. Cheminformatics_ , vol. 13, no. 1, p. 30, 2021.
  * 3. K. Huang, T. Fu, L. M. Glass, M. Zitnik, C. Xiao, and J. Sun, “DeepPurpose: a deep learning library for drug–target interaction prediction,” _Bioinformatics_ , vol. 36, no. 22-23, pp. 5545–5547, 2020.
  * 4. A. C. A. Nascimento, R. B. C. Prudêncio, and I. G. Costa, “A multiple kernel learning algorithm for drug-target interaction prediction,” _BMC Bioinformatics_ , vol. 17, no. 46, 2016.
  * 5. T. Nguyen, H. Le, T. P. Quinn, T. Nguyen, T. D. Le, and S. Venkatesh, “GraphDTA: predicting drug–target binding affinity with graph neural networks,” _Bioinformatics_ , vol. 37, no. 8, pp. 1140–1147, 2021.
  * 6. H. Öztürk, A. Özgür, and E. Ozkirimli, “DeepDTA: deep drug–target binding affinity prediction,” _Bioinformatics_ , vol. 34, no. 17, pp. i821–i829, 2018.
  * 7. S. Wang, P. Shan, Y. Zhao, and L. Zuo, “GanDTI: A multi-task neural network for drug-target interaction prediction,” _Comput. Biol. Chem._ , vol. 92, no. 9, p. 107476, 2021.
  * 8. Z. Wu, B. Ramsundar, E. N. Feinberg, J. Gomes, C. Geniesse, A. S. Pappu, K. Leswing, and V. Pande, “MoleculeNet: a benchmark for molecular machine learning,” _Chem. Sci._ , vol. 9, no. 2, pp. 513–530, 2018\.
  * 9. E. Vinogradova, A. Artykbayev, A. Amanatay, M. Karatayev, M. Mametkulov, A. Li, A. Suleimenov, A. Salimzhanov, K. Pats, R. Zhumagambetov, F. Molnár, V. Peshkov, and S. Fazli, “A biologically-inspired evaluation of molecular generative machine learning,” _arXiv_ , 2022.
  * 10. J. Tang, A. Szwajda, S. Shakyawar, T. Xu, P. Hintsanen, K. Wennerberg, and T. Aittokallio, “Making sense of large-scale kinase inhibitor bioactivity data sets: a comparative and integrative analysis,” _J. Chem. Inf. Model._ , vol. 54, no. 3, pp. 735–743, 2014.
  * 11. T. Liu, D. Tao, M. Song, and S. J. Maybank, “Algorithm-Dependent Generalization Bounds for Multi-Task Learning,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 39, no. 2, pp. 227–241, 2017.
  * 12. B. Ramsundar, S. Kearnes, P. Riley, D. Webster, D. Konerding, and V. Pande, “Massively Multitask Networks for Drug Discovery,” _arXiv_ , 2015\.
  * 13. Y. Zhang and Q. Yang, “A Survey on Multi-Task Learning,” _arXiv_ , 2021.
  * 14. Y. Zhang, “Multi-Task Learning and Algorithmic Stability,” _Proc. Conf. AAAI. Artif. Intell._ , vol. 29, no. 1, 2015.
  * 15. Z. Liu, X. Ye, X. Fang, F. Wang, H. Wu, and H. Wang, “Docking-based Virtual Screening with Multi-Task Learning,” _arXiv_ , 2021.
  * 16. S. Ruder, “An Overview of Multi-Task Learning in Deep Neural Networks,” Tech. Rep., 2017.
  * 17. G. E. Dahl, N. Jaitly, and R. Salakhutdinov, “Multi-task Neural Networks for QSAR Predictions,” Tech. Rep., 2014.
  * 18. J. Kirkpatrick, R. Pascanu, N. Rabinowitz, J. Veness, G. Desjardins, A. A. Rusu, K. Milan, J. Quan, T. Ramalho, A. Grabska-Barwinska, D. Hassabis, C. Clopath, D. Kumaran, and R. Hadsell, “Overcoming catastrophic forgetting in neural networks,” _Proc. Natl. Acad. Sci._ , vol. 114, no. 13, pp. 3521–3526, 2017.
  * 19. P. Mignone, G. Pio, S. Džeroski, M. Ceci, G. Pio, S. Džeroski, and M. Ceci, “Multi-task learning for the simultaneous reconstruction of the human and mouse gene regulatory networks,” _Sci. Rep._ , vol. 10, no. 1, p. 22295, 2020\.
  * 20. T. Liu, Y. Lin, X. Wen, R. N. Jorissen, and M. K. Gilson, “BindingDB: a web-accessible database of experimentally determined protein-ligand binding affinities,” _Nucleic Acid Res._ , vol. 35, pp. D198–D201, 2007.
  * 21. “Rdkit: Open-source cheminformatics software.” [Online]. Available: https://www.rdkit.org
  * 22. G. R. Bickerton, G. V. Paolini, J. Besnard, S. Muresan, and A. L. Hopkins, “Quantifying the chemical beauty of drugs,” _Nat. Chem._ , vol. 4, no. 2, pp. 90–98, 2012.
  * 23. M. Moret, L. Friedrich, F. Grisoni, D. Merk, and G. Schneider, “Generative molecular design in low data regimes,” _Nat. Mach. Intell._ , vol. 2, no. 3, pp. 171–180, 2020.
  * 24. M. Olivecrona, T. Blaschke, O. Engkvist, and H. Chen, “Molecular de-novo design through deep reinforcement learning,” _J. Cheminformatics_ , vol. 9, no. 1, p. 48, 2017.
  * 25. F. Grisoni, B. J. H. Huisman, A. L. Button, M. Moret, K. Atz, D. Merk, and G. Schneider, “Combining generative artificial intelligence and on-chip synthesis for de novo drug design,” _Sci. Adv._ , vol. 7, no. 24, p. eabg3338, 2021.