Generalized fractional calculus and some models of generalized counting
processes

Khrystyna Buchak

Faculty of Mathematics and Digital Technologies
Uzhgorod National University,
Narodna Square 3, 88000 Uzhgorod, Ukraine

Lyudmyla Sakhno

Faculty of Mechanics and Mathematics
Taras Shevchenko National University of Kyiv,
Volodymyrska St. 64, 01601 Kyiv, Ukraine

<EMAIL_ADDRESS>(K. Buchak),
<EMAIL_ADDRESS>(L. Sakhno)

###### Abstract

In the paper we consider models of generalized counting processes time-changed
by a general inverse subordinator, we characterize their distributions and
present governing equations for them. The equations are given in terms of the
generalized fractional derivatives, namely, convolutions-type derivatives with
respect to Bernštein functions. Some particular examples are presented.

Key words: time-change, Poisson process, generalized counting process,
subordinator, inverse subordinator, generalized fractional derivatives

2010 MSC: 60G50; 60G51; 60G55

##  1 Introduction

Stochastic processes with random time-change has become a well-established and
highly ramified branch of the modern theory of stochastic processes which
attracts more and more attention due to various applications in financial,
biological, ecological, physical, technical and other fields of research.

Rich class of models is provided by time-changed Poisson processes, of which
the most intensively studied are two fractional extensions of the Poisson
process, namely, the space-fractional and time-fractional Poisson processes
obtained by choosing a stable subordinator or its inverse in the role of time
correspondingly. We refer, for example, to papers [8], [18], [20], [1], [2],
[16], [19], among many others (see also references therein). In particular, in
paper [20] a general class of time-changed Poisson processes
$N^{f}(t)=N(H^{f}(t))$, $t>0$, was introduced and studied, where $N(t)$ is a
Poisson process and $H^{f}(t)$ is an arbitrary subordinator with Laplace
exponent $f,$ independent of $N(t)$, their distributional properties, hitting
times and governing equations were presented (see, also [8]). In paper [5]
Poisson processes time changed by general inverse subordinators were studied,
the governing equations for their marginal distributions were presented and
some other properties were described. Poisson process itself, being in a sense
a core object concerning applicability to count data and simple tractability,
however, as a reverse side of its simplicity, is a rather restrictive model.
Therefore, it has been quite natural to search for its extensions and
generalizations to provide some new useful features and properties needed for
applications. For example, time changed processes $N(H^{f}(t))$ allow for
jumps of arbitrary size and other interesting properties ([20], [8]).

Recently, generalized counting processes, their time-changed versions and
fractional extensions have been intensively investigated, including Poisson
and Skellam processes of order $k$, Pólya-Aeppli process of order $k$, Bell-
Touchard process, Poisson-logarithmic process, etc. and their fractional
extensions. We refer to papers [7], [15], [6], [10], [11], [12], [13], to
mention only few, see also references therein.

In the present paper we consider several particular models of generalized
counting processes time-changed by a general inverse subordinator,
characterize their distributions and present governing equations for them,
which are obtained following the technique presented in [5]. The equations are
given in terms of the generalized Caputo-Djrbashian derivatives, which called
also convolutions-type derivatives with respect to Bernštein functions. These
generalized derivatives were introduced in [14] and [21] and have been widely
used to describe various advanced stochastic models.

The convolution-type derivatives allow to study the properties of
subordinators and their inverses in the unifying manner ([21, 17]), in
particular, the governing equations for their densities can be given in terms
of convolution-type derivatives. These properties will be the key tools for
our study: by considering several models of time-changed processes, we
elucidate the technique leading from the equations for the densities of
inverse subordinators and their Laplace transforms to the equations for
probabilities of processes time-changed by inverse subordinators and equations
for some related functions.

The paper is organized as follows. Section 2 collects the basic definitions
and facts on the convolution-type derivatives needed in our study. We also
consider, as a warmup example, the time changed Poisson process
$N(Y^{f}(t))$,$t>0$, where $Y^{f}$ is the inverse to the subordinator with
Laplace exponent $f,$ independent of $N(t)$. We partly extend the results from
[5] for this process. This simplest case allows to demonstrate very
transparently the technique underpinning the derivations of the corresponding
results for more general models considered in the next sections. Namely, in
Sections 3 and 4 we study the processes
$N^{\psi,f}(t)=N\left(H^{\psi}(Y^{f}(t))\right)$, $t>0$, and
$M^{\psi,f}(t)=M\left(H^{\psi}(Y^{f}(t))\right)$, $t>0$, where $M$ is the
generalized counting process. Time change is done by the independent
subordinator $H^{\psi}$ and inverse subordinator $Y^{f}$ related to Bernštein
functions $\psi$ and $f$ correspondingly. We present the differential-
difference equations for probabilities of the processes $N^{\psi,f}$ and
$M^{\psi,f}$, which are given in terms of the convolutions-type derivative
with respect to function $f$ and the difference operator related to function
$\psi$. Expressions for the probabilities of these processes are also
presented in terms of the Laplace transform of the process $Y^{f}$. Some other
properties of these processes are described and a comparison with related
results in the literature are given. As examples, we consider the processes
$N^{\psi,f}$ and $M^{\psi,f}$, where $H^{\psi}$ is a compound Poisson-Gamma
process.

##  2 Preliminaries

Generalized fractional derivatives. We review the main definitions and some
facts on the generalized fractional derivatives, which will be important for
our study. (For more details see, e.g., [17, 21].)

Let $f(x)$ be a Bernštein function:

$f(x)=a+bx+\int_{0}^{\infty}\left(1-e^{-xs}\right)\overline{\nu}(ds),\,\,x>0,\,\,a,b\geq
0,$ (1)

$\overline{\nu}(ds)$ is a non-negative measure on $(0,\infty)$ (the Lévy
measure for $f(x)$) such that

$\int_{0}^{\infty}\left(s\wedge 1\right)\overline{\nu}(ds)<\infty.$

The generalized Caputo-Djrbashian (C-D) derivative, or convolution-type
derivative, with respect to the Bernštein function $f$ is defined on the space
of absolutely continuous functions as follows ([21], Definition 2.4):

$\mathcal{D}^{f}_{t}u(t)=b\frac{d}{dt}u(t)+\int_{0}^{t}\frac{\partial}{\partial
t}u(t-s)\nu(s)ds,$ (2)

where $\nu(s)=a+\overline{\nu}(s,\infty)$ is the tail of the Lévy measure
$\overline{\nu}(s)$ of the function $f$.

In the case where $f(x)=x^{\alpha},x>0,\alpha\in(0,1)$, the derivative (2)
reduces to the fractional C-D derivative:

$\mathcal{D}^{f}_{t}u(t)=\frac{d^{\alpha}}{dt^{\alpha}}u(t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}\frac{u^{\prime}(s)}{\left(t-s\right)^{\alpha}}ds.$
(3)

For the Laplace transform of the derivative (2) the following relation holds
([21], Lemma 2.5):

$\mathcal{L}\left[\mathcal{D}^{f}_{t}u\right](s)=f(s)\mathcal{L}\left[u\right](s)-\frac{f(s)}{s}u(0),s>s_{0},$

for $u$ such that $|u(t)|\leq\mathsf{M}e^{s_{0}t}$, $M$ and $s_{0}$ are some
constants. Similarly to the C-D fractional derivative, the convolution type
derivative can be alternatively defined via its Laplace transform.

The generalization of the classical Riemann-Liouville (R-L) fractional
derivative is introduced in [21] by means of another convolution-type
derivative with respect to $f$ given by the following formula:

$\mathbb{D}_{t}^{f}u(t)=b\frac{d}{dt}u(t)+\frac{d}{dt}\int_{0}^{t}u(t-s)\nu(s)ds.$
(4)

The derivatives $\mathcal{D}^{f}_{t}$ and $\mathbb{D}_{t}^{f}$ are related as
follows (see, [21], Proposition 2.7):

$\mathbb{D}_{t}^{f}u(t)=\mathcal{D}^{f}_{t}u(t)+\nu(t)u(0).$

Bernštein functions are associated in a natural way with subordinators.

Let $H(t),t\geq 0,$ be a subordinator, that is, nondecreasing Lévy process.
Its Laplace transform is of the form:

$\mathcal{L}[H(t)](s)=\mathsf{E}e^{-sH(t)}=e^{-tf(s)},$

where the function $f$, called the Laplace exponent, is a Bernštein function.
Consider a subordinator $H^{f}$ with the Laplace exponent $f$ given by (1),
and let $Y^{f}$ be its inverse process defined as

$Y^{f}(t)=\inf\left\\{s\geq 0:H^{f}(s)>t\right\\}.$ (5)

It was shown in [21] that the distribution of the inverse process $Y^{f}$ has
a density $\ell_{f}(t,x)={\rm P}\\{Y^{f}(t)\in dx\\}/dx$ provided that the
following condition holds:

Condition I. $\overline{\nu}(0,\infty)=\infty$ and the tail
$\nu(s)=a+\overline{\nu}(s,\infty)$ is absolutely continuous.

The Laplace transform of the density with respect to $t$ has the form ([21]):

$\mathcal{L}_{t}\left(\ell_{f}(t,x)\right)(r)=\frac{f(r)}{r}e^{-xf(r)}.$

The density $\ell_{f}(t,u)$ of the inverse process $Y^{f}$ satisfies the
following equation ([21], Theorem 4.1):

${\mathbb{D}}_{t}^{f}\ell_{f}(t,u)=-\frac{\partial}{\partial u}\ell_{f}(t,u),$
(6)

subject to

$\ell_{f}(t,{u}/{b})=0,\,\,\ell_{f}(t,0)=\nu(t),\,\,\ell_{f}(0,u)=\delta(u).$
(7)

The space Laplace transform of the density $\ell_{f}(t,x)$

$\tilde{\ell_{f}}(t,\lambda)=\int_{0}^{\infty}e^{-\lambda
x}\ell_{f}(t,x)dx={\mathsf{E}}e^{-\lambda L^{f}(t)}$ (8)

is an eigenfunction of the operator $\mathcal{D}^{f}_{t}$, that is, satisfies
the equation

$\mathcal{D}^{f}_{t}\tilde{\ell_{f}}(t,\lambda)=-\lambda\tilde{\ell_{f}}(t,\lambda).$
(9)

(see, [5, 14, 17]). In the case where $f(x)=x^{\alpha},x>0,\alpha\in(0,1)$,
$\tilde{\ell_{f}}(t,\lambda)=\mathcal{E}_{\alpha}(-\lambda t^{\alpha})$, where
$\mathcal{E}_{\alpha}(\cdot)$ is the Mittag-Leffler function.

Equations (6)–(9) are important for the study of the processes time changed by
inverse subordinators as can be seen in what follows. We first consider the
Poisson process with an inverse subordinator as a simplest example which
demonstrates the technique to be applied further to more general models.

Poisson process time changed by an inverse subordinator. Let $N(t)$ be the
Poisson process with intensity $\lambda,$ and $Y^{f}(t)$ be the inverse
subordinator.

######  Theorem 1.

Let Condition I hold. Then the marginal distributions
$p_{k}^{f}(t)=P\left\\{N(Y^{f}(t))=k\right\\}$, $k=0,1,\dots$ satisfy the
differential equation

$\mathcal{D}^{f}_{t}p_{k}^{f}(t)=\lambda\left(p_{k}^{f}(t)-p_{k-1}^{f}(t)\right)$
(10)

and

$p_{k}^{f}(t)=\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}\tilde{\ell}_{f}(t,\lambda),$
(11)

where $\tilde{\ell}_{f}(t,\lambda)=\int_{0}^{\infty}e^{-\lambda
u}\ell_{f}(t,u)du$ is the Laplace transform of the density of the inverse
subordinator $Y^{f}(t).$

The probability generating function of the process $N^{f}$ has the form

$G^{f}(u,t)=\tilde{\ell}(t,\lambda(1-u)),\,|u|<1,$ (12)

and satisfies the equation

$\mathcal{D}^{f}_{t}G^{f}(u,t)=-\lambda(1-u)G^{f}(u,t)$ (13)

with $G^{f}(u,0)=1$.

###### Proof.

Proof of equation (10). For the probabilities $p_{x}^{f}(t)$ we have:

$\displaystyle p_{x}^{f}(t)$ $\displaystyle=$
$\displaystyle\mathsf{P}\left\\{N\left(Y^{f}(t)\right)=x\right\\}=\int_{0}^{\infty}p_{x}(u)\ell_{f}(t,u)du$
(14) $\displaystyle=$ $\displaystyle\int_{0}^{\infty}\frac{e^{-\lambda
u}(\lambda u)^{x}}{x!}\ell_{f}(t,u)du,\quad x=0,1,2,\dots$

We take the generalized R-L convolution-type derivative $\mathbb{D}^{f}_{t}$
given by (4) and use the equation for the density $\ell^{f}(t,u)$ of the
inverse subordinator:

${\mathbb{D}}^{f}_{t}\ell_{f}(t,u)=-\frac{\partial}{\partial u}\ell_{f}(t,u),$
(15)

with

$\ell_{f}(t,0)=\nu(t),\,\,\,\,\ell_{f}(0,u)=\delta(u).$ (16)

We obtain:

$\displaystyle{\mathbb{D}}^{f}_{t}p_{x}^{f}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}p_{x}(u){\mathbb{D}}^{f}_{t}\ell_{f}(t,u)du=-\int_{0}^{\infty}p_{x}(u)\frac{\partial}{\partial
u}\ell_{f}(t,u)du$ (17) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\ell_{f}(t,u)\frac{\partial}{\partial
u}p_{x}(u)du-p_{x}(u)\ell_{f}(t,u)\big{|}_{u=0}^{\infty}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\ell_{f}(t,u)(-\lambda[p_{x}(u)-p_{x-1}(u)])du+p_{x}(0)\ell_{f}(t,0)$
$\displaystyle=$
$\displaystyle-\lambda\left[p^{f}_{x}(t)-p^{f}_{x-1}(t)\right]+p_{x}(0)\nu(t).$

Using the relation

$\mathbb{D}^{f}_{t}u(t)=\mathcal{D}^{f}_{t}+\nu(t)u(0)$ (18)

between convolution derivatives of C-D and R-L types, we have:

$\mathcal{D}^{f}_{t}p_{x}^{f}(t)=\mathbb{D}^{f}_{t}p_{x}^{f}(t)-\nu(t)p_{x}^{f}(0),$
(19)

note also that

$p_{x}^{f}(0)=\int_{0}^{\infty}p_{x}(u)\ell_{f}(0,u)du=\int_{0}^{\infty}p_{x}(u)\delta(u)du=p_{x}(0)=1.$
(20)

From (17), taking into account (19)-(20), we finally obtain:

$\mathcal{D}^{f}_{t}p_{x}^{f}(t)=-\lambda\left[p_{x}^{f}(t)-p_{x-1}^{f}(t)\right].$

Proof of equation (11):

$\displaystyle\mathsf{P}\left\\{N(Y^{f}(t))=k\right\\}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathsf{P}(N(s)=k)\mathsf{P}(Y^{f}(t)\in ds)$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}e^{-\lambda
s}\mathsf{P}\left\\{Y^{f}(t)\in ds\right\\}$ $\displaystyle=$
$\displaystyle\mathsf{E}\left[\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}\exp\left\\{-\lambda
Y^{f}(t)\right\\}\right]$ $\displaystyle=$
$\displaystyle\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}\mathsf{E}\exp\left(-\lambda
Y^{f}(t)\right)=\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}\tilde{\ell}_{f}(t,\lambda).$

We used above that the distribution of the Poisson process $N$ with the rate
$\lambda$ can be writen as follows:

$\displaystyle
p_{k}(t)=\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}e^{-\lambda
t}.$

For the probability generating function we have

$G^{f}(u,t)=\mathsf{E}u^{N(Y^{f}(t))}\int_{0}^{\infty}e^{-s\lambda(1-u)}\ell_{f}(t,s)ds=\mathsf{E}e^{-\lambda(1-u)Y^{f}(t)}=\tilde{\ell}_{f}(t,\lambda(1-u)),$

therefore, we obtain formula (12). Now we use the fact, that
$\tilde{\ell}_{f}(t,\lambda)$ is an eigenfunction of the operator
$\mathcal{D}^{f}_{t}$ (see (9)), from which we conclude that $G^{f}(u,t)$,
being given by $\tilde{\ell}_{f}(t,\lambda(1-u))$, satisfies equation (13). ∎

######  Remark 1.

Note that (9) can be deduced directly from equation (10) by taking $k=1$,
since
$p_{0}^{f}(t)=\mathsf{P}\left\\{N\left(Y^{f}(t)\right)=0\right\\}=\int_{0}^{\infty}e^{-\lambda
u}\ell_{f}(t,u)du=\tilde{\ell}(t,\lambda)$.

##  3 Models of time-changed Poisson processes

Consider the time-changed Poisson processes
$N^{\psi}(t)=N\left(H^{\psi}(t)\right)$, $t>0$, where $N(t)$ is the Poisson
process with intensity $\lambda$ and $H^{\psi}(t)$ is the subordinator with
Bernštein function $\psi(u)$, independent of $N(t)$. This class of processes
was introduced and studied in [20] and called by the authors Poisson processes
with Bernštein intertimes. It was shown in [20] that the distributions of
$N^{\psi}(t)$, $t>0$, can be presented as follows:

$p_{k}^{\psi}(t)=\mathsf{P}\left\\{N^{\psi}(t)=k\right\\}=\frac{(-1)^{k}}{k!}\frac{d^{k}}{dt^{k}}e^{-t\psi(\lambda
u)}\arrowvert_{u=1},$ (21)

and satisfy the difference-differential equations:

$\frac{d}{dt}p_{k}^{\psi}(t)=-\psi(\lambda)p_{k}^{\psi}(t)+\sum_{m=1}^{k}\frac{\lambda^{m}}{m!}p_{k-m}^{\psi}(t)\int_{0}^{\infty}e^{-s\lambda}s^{m}\nu(ds),\,k\geq
0,\,t>0,$

which can be also written as

$\frac{d}{dt}p_{k}^{\psi}(t)=-\psi(\lambda)p_{k}^{\psi}(t)-\sum_{m=1}^{k}\psi^{(m)}(\lambda)\frac{(-\lambda)^{m}}{m!}p_{k-m}^{\psi}(t),\,k\geq
0,\,t>0,$ (22)

with the usual initial conditions: $p_{0}^{\psi}(0)=1$, $p_{k}^{\psi}(0)=0$
for $k\geq 1$. The last equation can be represented in the form:

$\frac{d}{dt}p_{k}^{\psi}(t)=-\psi(\lambda(I-B))p_{k}^{\psi}(t),\,k\geq
0,\,t>0,$ (23)

where $B$ is the backshift operator: $Bp_{k}^{f}(t)=p_{k-1}^{f}(t)$, and it is
supposed that $p_{-1}(t)=0$.

The probability generating function of the process $N^{\psi}$ is of the form
([20]):

$G^{\psi}(u,t)=\mathsf{E}u^{N^{\psi}(t)}=e^{-t\psi(\lambda(1-u))},\,|u|<1.$
(24)

We refer for more detail on these processes to [20].

Consider the process $N^{\psi}$ time-changed by an independent inverse
subordinator $Y^{f}$. We can state the following result.

######  Theorem 2.

The process $N^{\psi,f}(t)=N\left(H^{\psi}(Y^{f}(t))\right)$, $t>0$, has
probability distribution function

$p_{k}^{\psi,f}(t)=\mathsf{P}\left\\{N^{\psi,f}(t)=k\right\\}=\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}\tilde{\ell}_{f}(t,\psi(\lambda)),$
(25)

and the probabilities $p_{k}^{\psi,f}$ satisfy the following equation

$\mathcal{D}^{f}_{t}p_{k}^{\psi,f}(t)=-\psi\left(\lambda\left(I-B\right)\right)p_{k}^{\psi,f}(t)$
(26)

with the usual initial conditions.

The probability generating function of the process $N^{\psi,f}$ has the form

$G^{\psi,f}(u,t)=\tilde{\ell}(t,\psi(\lambda(1-u))),\,|u|<1,$ (27)

and satisfies the equation

$\mathcal{D}^{f}_{t}G^{\psi,f}(u,t)=-\psi(\lambda(1-u))G^{\psi,f}(u,t)$ (28)

with $G^{\psi,f}(u,0)=1$.

###### Proof.

To prove equation (25) we perform the calculations as follows:

$\displaystyle p_{k}^{\psi,f}(t)$ $\displaystyle=$
$\displaystyle\mathsf{P}\left\\{N\left(H^{\psi}(Y^{f}(t))\right)=k\right\\}$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathsf{P}\left\\{N(s)=k\right\\}\mathsf{P}\left\\{H^{\psi}(Y^{f}(t))\in
ds\right\\}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}e^{-\lambda
s}\mathsf{P}\left\\{H^{\psi}(Y^{f}(t))\in ds\right\\}$ $\displaystyle=$
$\displaystyle\mathsf{E}\left[\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}\exp\left\\{-\lambda
H^{\psi}(Y^{f}(t))\right\\}\right]$ $\displaystyle=$
$\displaystyle\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}\mathsf{E}\exp\left\\{-\lambda
H^{\psi}(Y^{f}(t))\right\\}.$

Now it is left to note that

$\mathsf{E}\exp\left\\{-\lambda
H^{\psi}(Y^{f}(t))\right\\}=\mathsf{E}\exp\left(-\psi(\lambda)Y^{f}(t)\right)=\tilde{\ell}_{f}(t,\psi(\lambda)).$

Let us state equation (26). We have

$p_{k}^{\psi,f}(t)=\mathsf{P}\left\\{N\left(H^{\psi}(Y^{f}(t))\right)=k\right\\}=\int_{0}^{\infty}p_{k}^{\psi}(u)\ell_{f}(t,u)du,\quad
k=0,1,2,\dots$

We repeat the same lines as in the proof of Theorem 1. We take the generalized
R-L derivative $\mathbb{D}_{t}^{f}$ and use equations (22)–(23):

$\displaystyle\mathbb{D}_{t}^{f}p_{k}^{\psi,f}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}p_{k}^{\psi}(u)\mathbb{D}_{t}^{f}\ell_{f}(t,u)du=-\int_{0}^{\infty}p_{k}^{\psi}(u)\frac{d}{du}\ell_{f}(t,u)du$
(29) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\ell_{f}(t,u)\frac{d}{du}p_{k}^{\psi}(u)du-
p_{k}^{\psi}(u)\ell_{f}(t,u)\big{|}_{0}^{\infty}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\ell_{f}(t,u)\left(-\psi(\lambda(I-B))\right)p_{k}^{\psi}(u)du+p_{k}^{\psi}(0)\ell_{f}(t,0)$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\ell_{f}(t,u)\Big{(}-\Big{(}\psi(\lambda)p_{k}^{\psi}(u)$
$\displaystyle+\sum_{m=1}^{k}\psi^{(m)}(\lambda)\frac{(-\lambda)^{m}}{m!}p_{k-m}^{\psi}(u)\Big{)}\Big{)}du+p_{k}^{\psi}(0)\nu(t)$
$\displaystyle=$
$\displaystyle-\psi(\lambda)p_{k}^{\psi,f}(t)-\sum_{m=1}^{k}\psi^{(m)}(\lambda)\frac{(-\lambda)^{m}}{m!}p_{k-m}^{\psi,f}(t)+p_{k}^{\psi}(0)\nu(t)$
$\displaystyle=$
$\displaystyle-\psi(\lambda(I-B))p_{k}^{\psi,f}(t)+p_{k}^{\psi}(0)\nu(t).$

Using the relation between generalized derivatives of C-D and R-L types, we
can write

$\mathcal{D}^{f}_{t}p_{k}^{\psi,f}(t)=\mathbb{D}_{t}^{f}p_{k}^{\psi,f}(t)-\nu(t)p_{k}^{\psi,f}(0);$

and we have (see (7)):

$p_{k}^{\psi,f}(0)=\int_{0}^{\infty}p_{k}^{\psi}(u)\ell_{f}(0,u)du=\int_{0}^{\infty}p_{k}^{\psi}(u)\delta(u)du=p_{k}^{\psi}(0)=1.$

Therefore, we obtain

$\mathcal{D}^{f}_{t}p_{k}^{\psi,f}(t)=-\psi(\lambda(I-B))p_{k}^{\psi,f}(t).$

Using the expression for $G^{\psi}(u,t)$ given by (24), we calculate
$G^{\psi,f}:$

$\displaystyle G^{\psi,f}(u,t)$ $\displaystyle=$
$\displaystyle\mathsf{E}u^{N\left(H^{\psi}(Y^{f}(t))\right)}=\int_{0}^{\infty}G^{\psi}(u,s)\ell_{f}(t,s)ds=$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}e^{-s\psi(\lambda(1-u))}\ell_{f}(t,s)ds=\mathsf{E}e^{-\psi(\lambda(1-u))Y^{f}(t)}=$
$\displaystyle=$ $\displaystyle\tilde{\ell}_{f}(t,\psi(\lambda(1-u))),$

and since $\tilde{\ell}_{f}(t,\psi(\lambda(1-u)))$ is an eigenfunction of
$\mathcal{D}^{f}_{t}$ with the eigenvalue $\psi(\lambda(1-u))$ (see (9)), it
follows that $G^{\psi,f}(u,t)$ satisfies equation (28). ∎

######  Remark 2.

We notice at once that equations (26) for the probabilities of the process
$N^{\psi,f}(t)=N(H^{\psi}(Y^{f}(t)))$ mimic the corresponding equations for
the process $N^{\psi}(t)=N(H^{\psi}(t))$, that is, the process before the time
change by an inverse subordinator, only the ordinary derivative in time is
changed for the generalized fractional derivative. This can be anticipated
quite straightforwardly, in view of the technique used for the proof, and
holds also for other models, like those below and in the next section. On the
other side, the equations (26) tell us, that for the probabilities of the
processes with double time change like $N\left(H^{\psi}(Y^{f}(t))\right)$, the
action (in time) of the operator $\mathcal{D}^{f}_{t}$ (which is related to
$Y^{f}$) is equal to the action (in space) of the operator
$-\psi\left(\lambda\left(I-B\right)\right)$ (which is related to $H^{\psi}$,
and depends also of the outer process $N$). Some more insight on these
operators can be found within the approach applied in the paper [1]. Let
$Y^{f}=Y^{\beta}$ be the inverse stable subordinator with a parameter $\beta$,
that is, $f(\lambda)=\lambda^{\beta}$, then equation (26) becomes

$\mathcal{D}^{\beta}_{t}p_{k}^{\psi,\beta}(t)=-\psi\left(\lambda\left(I-B\right)\right)p_{k}^{\psi,\beta}(t),$
(30)

where $\mathcal{D}^{\beta}_{t}$ is the Caputo-Djrbashian fractional derivative
(3). If we suppose furthermore that $H^{\psi}=H^{\alpha}$ is the stable
subordinator with a parameter $\alpha$, that is,
$\psi(\lambda)=\lambda^{\alpha}$, then we come to the equation

$\mathcal{D}^{\beta}_{t}p_{k}^{\alpha,\beta}(t)=-\lambda^{\alpha}\left(I-B\right)^{\alpha}p_{k}^{\alpha,\beta}(t).$
(31)

Equations (30) and (31) were stated in paper [1], where the particular
representation of the operator $\psi\left(\lambda\left(I-B\right)\right)$ was
used, and the equations for the probabilities of time changed processes were
stated using the interplay between the operator
$\psi\left(\lambda\left(I-B\right)\right)$ and the fractional operator
$\mathcal{D}^{\beta}_{t}$ in the Fourier domain. The approach in [1] can be
extended for the case where we deal with the generalized fractional operator
$\mathcal{D}^{f}_{t}$, as soon as we know that its eigenfunction is given by
$\tilde{\ell}_{f}(t,\lambda)$ (see (9)). We refer for more detail on this
approach to [1].

Consider a generalization of counting processes with Bernštein intertimes,
which was introduced in the paper [8]:

$N^{\psi_{1},\psi_{2},\dots,\psi_{n}}(t)=N\left(\sum_{j=1}^{n}H^{\psi_{j}}(t)\right),t\geq
0,$

where $H^{\psi_{i}},i=1,\dots n$ are $n$ independent subordinators with
Bernštein functions $\psi_{i}$, independent of the Poisson process $N$. In [8]
the following result was stated.

######  Proposition 1 ([8]).

The distribution of the subordinated process
$N^{\psi_{1},\psi_{2},\dots,\psi_{n}}(t),$ $t\geq 0$, is the solution to the
Cauchy problem

$\frac{d}{dt}p_{k}^{\psi_{1},\psi_{2},\dots,\psi_{n}}(t)=-\sum_{j=1}^{n}\psi_{j}(\lambda(I-B))p_{k}^{\psi_{1},\psi_{2},\dots,\psi_{n}}(t),$

with the usual initial conditions

$p_{k}^{\psi_{1},\psi_{2},\dots,\psi_{n}}(0)=\begin{cases}1,\quad k>0,\\\
0,\quad k=0.\\\ \end{cases}$

Consider the process
$N^{\psi_{1},\psi_{2},\dots,\psi_{n},f}(t)=N^{\psi_{1},\psi_{2},\dots,\psi_{n}}(Y^{f}(t)),t\geq
0.$

######  Theorem 3.

The process $N^{\psi_{1},\psi_{2},\dots,\psi_{n},f}$ has probability
distribution function

$p_{k}^{\psi_{1},\psi_{2},\dots,\psi_{n},f}(t)=\frac{\left(-\lambda\partial_{\lambda}\right)^{k}}{k!}\tilde{\ell}_{f}(t,\sum_{j=1}^{n}\psi_{j}(\lambda)),$

and the probabilities $p_{k}^{\psi_{1},\psi_{2},\dots,\psi_{n},f}(t)$ satisfy
the following equation

$\mathcal{D}^{f}_{t}p_{k}^{\psi_{1},\psi_{2},\dots,\psi_{n},f}(t)=-\sum_{j=1}^{n}\psi_{j}\left(\lambda\left(I-B\right)\right)p_{k}^{\psi_{1},\psi_{2},\dots,\psi_{n},f}(t)$

with the usual initial conditions.

Proof of Theorem 3 is obtained by the same reasoning as that of Theorem 2,
using Proposition 1.

Example 1. Consider the time-changed process $N^{GN}(t)=N_{1}(G_{N}(t)),$
$t>0$, where $N_{1}(t)$ is the Poisson process with intensity $\lambda_{1},$
and $G_{N}(t)=G(N(t))$, $t>0$, is the compound Poisson-Gamma subordinator with
parameters $\lambda,\alpha,\beta$, that is, with the Laplace exponent
$\psi_{GN}(u)=\lambda\beta^{\alpha}\left(\beta^{-\alpha}-(\beta+u)^{-\alpha}\right).$
In the case when $\alpha=1$ we have the compound Poisson-exponential
subordinator, which we will denote as $E_{N}(t)$, and the corresponding time-
changed process as $N^{E}(t)$. The detailed study of the processes $N^{GN}(t)$
and $N^{E}(t)$ was presented in [3], [4].

We now consider the processes $N^{GN,f}(t)=N^{GN}(Y^{f}(t))$ and
$N^{E,f}(t)=N^{E}(Y^{f}(t))$, $t\geq 0.$

The next theorem is obtained as a corollary of Theorem 2.

######  Theorem 4.

The process $N^{GN,f}(t)$ has probability distribution function

$p_{k}^{GN,f}(t)=\frac{\left(-\lambda_{1}\partial_{\lambda_{1}}\right)^{k}}{k!}\tilde{\ell}_{f}(t,\psi_{GN}(\lambda_{1})),$

and $p_{k}^{GN,f}(t)$ satisfy the following equation

$\displaystyle\mathcal{D}^{f}_{t}p_{k}^{GN,f}(t)$ $\displaystyle=$
$\displaystyle\left(\frac{\lambda\beta^{\alpha}}{(\lambda_{1}+\beta)^{\alpha}}-\lambda\right)p_{k}^{GN,f}(t)$
$\displaystyle+\frac{\lambda\beta^{\alpha}}{(\lambda_{1}+\beta)^{\alpha}}\sum_{m=1}^{k}\frac{\lambda_{1}^{m}}{\left(\lambda_{1}+\beta\right)^{m}}\frac{\Gamma(m+\alpha)}{m!\Gamma(\alpha)}p_{k-m}^{GN,f}(t).$

The process $N^{E,f}(t)=N_{1}(E_{N}(Y^{f}(t)))$ has probability distribution
function

$p_{k}^{E,f}(t)=\frac{\left(-\lambda_{1}\partial_{\lambda_{1}}\right)^{k}}{k!}\tilde{\ell}_{f}(t,\psi_{E}(\lambda_{1})),$

where $\psi_{E}(u)=\frac{\lambda u}{\beta+u},$ and $p_{k}^{E,f}(t)$ satisfy
the following equation

$\mathcal{D}^{f}_{t}p_{k}^{E,f}(t)=-\lambda\frac{\lambda_{1}}{\lambda_{1}+\beta}p_{k}^{E,f}(t)+\frac{\lambda\beta}{\lambda_{1}+\beta}\sum_{m=1}^{k}\left(\frac{\lambda_{1}}{\lambda_{1}+\beta}\right)^{m}p_{k-m}^{E,f}(t).$

##  4 Models of time-changed generalized counting processes

Consider the generalized counting process (GCP) $M(t)$, $t\geq 0$, which was
introduced in [7]. The probabilities
$\tilde{p}_{n}(t)=\mathsf{P}\left\\{M(t)=n\right\\}$ depend on $k$ parameters
$\lambda_{1},\ldots,\lambda_{k}$ and are given by the formula

$\tilde{p}_{n}(t)=\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\left(\lambda_{j}t\right)^{x_{j}}}{x_{j}!}e^{-\Lambda
t},n\geq 0,$ (32)

where
$\Omega(k,n)=\left\\{(x_{1},\dots,x_{k})\,:\,\sum_{j=1}^{k}jx_{j}=n,x_{j}\in
N_{0}\right\\}$, $\Lambda=\sum_{j=1}^{k}\lambda_{j}$.
GCP performs $k$ kinds of jumps of amplitude $1,2,\ldots,k$ with rates
$\lambda_{1},\ldots,\lambda_{k}$.

The probability generating function of GCP is given by ([11]):

$\tilde{G}(u,t)=\mathsf{E}u^{M(t)}=\exp\Big{\\{}-\sum_{j=1}^{k}\lambda_{j}(1-u^{j})t\Big{\\}},\,|u|<1.$
(33)

The probabilities (32) satisfy

$\frac{\tilde{p}_{n}(t)}{dt}=-\Lambda\tilde{p}_{n}(t)+\sum_{j=1}^{\min\\{n,k\\}}\lambda_{j}\tilde{p}_{n-k}(t),\,n\geq
0,$

with the usual initial condition. The probabilities $\tilde{p}_{n}(t)$ can be
also written as

$\tilde{p}_{n}(t)=\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\left(-\partial_{\Lambda}\right)^{z_{k}}e^{-\Lambda
t},\,n\geq 0,$ (34)

where $z_{k}=\sum_{j=1}^{k}x_{j}$ (see [13]).

Recently various models of time-changed GCP and fractional GCP were studied,
in particular, with time-change given by Lévy subordinators, inverse
subordinators, including the cases of some specific subordinators, and also
the fractional extensions of processes of the form $M(H^{\psi}(t))$ (see, for
example, [7], [10], [11], [12], [13], and references therein).

Following the lines of the previos section, we consider the time-changed
process $\mathsf{M}^{\psi,f}(t)=M(H^{\psi}(Y^{f}(t))$, that is, with double
time-change by an independent subordinator $H^{\psi}$ and an inverse
subordinator $Y^{f}$, which are independent of $M$. To the best of our
knowledge, such general case has not been presented in the literature. In the
next theorem we characterize its probabilities
$\tilde{p}_{n}^{\psi,f}(t)=\mathsf{P}\left\\{\mathsf{M}^{\psi,f}(t)=n\right\\}$
and probability generating function.

######  Theorem 5.

The process $\mathsf{M}^{\psi,f}$ has probability distribution function

$\tilde{p}_{n}^{\psi,f}(t)=\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\left(-\partial_{\Lambda}\right)^{z_{k}}\tilde{\ell}_{f}(t,\psi(\Lambda)),\,n\geq
0,$ (35)

and $\tilde{p}_{n}^{\psi,f}(t)$ satisfy the following equation

$\displaystyle\mathcal{D}^{f}_{t}\tilde{p}_{n}^{\psi,f}(t)=-\psi(\Lambda)\tilde{p}_{n}^{\psi,f}(t)-\sum_{m=1}^{n}\sum_{\Omega(k,m)}\psi^{(z_{k})}\left(\Lambda\right)\prod_{j=1}^{k}\frac{(-\lambda_{j})^{x_{j}}}{x_{j}!}\tilde{p}_{n-m}^{\psi,f}(t),$
(36)

with initial conditions

$\tilde{p}_{n}^{\psi,f}(0)=\begin{cases}1,\quad n=0,\\\ 0,\quad n\geq 1.\\\
\end{cases}$

The probability generating function of the process $\mathsf{M}^{\psi,f}$ is of
the form

$\tilde{G}^{\psi,f}(u,t)=\tilde{\ell}_{f}\Big{(}t,\psi\Big{(}\sum_{j=1}^{k}\lambda_{j}(1-u^{j})\Big{)}\Big{)},\,|u|<1,$
(37)

and satisfies the equation

$\mathcal{D}^{f}_{t}\tilde{G}^{\psi,f}(u,t)=-\psi\Big{(}\sum_{j=1}^{k}\lambda_{j}(1-u^{j})\Big{)}\tilde{G}^{\psi,f}(u,t)$
(38)

with $\tilde{G}^{\psi,f}(u,0)=1$.

###### Proof.

For calculating $\tilde{p}_{n}^{\psi,f}(t)$ we use the formula (34):

$\displaystyle\tilde{p}_{n}^{\psi,f}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\tilde{p}_{n}(s)P\left\\{H^{\psi}(Y^{f}(t))\in
ds\right\\}=$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\left(-\partial_{\Lambda}\right)^{z_{k}}e^{-\Lambda
t}P\left\\{H^{\psi}(Y^{f}(t))\in ds\right\\}=$ $\displaystyle=$
$\displaystyle\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\mathsf{E}\left(-\partial_{\Lambda}\right)^{z_{k}}\exp\left\\{-\Lambda
H^{\psi}(Y^{f}(t))\right\\}=$ $\displaystyle=$
$\displaystyle\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\left(-\partial_{\Lambda}\right)^{z_{k}}\mathsf{E}\exp\left\\{-\Lambda
H^{\psi}(Y^{f}(t))\right\\}.$

This implies formula (35), since

$\mathsf{E}\exp\left\\{-\Lambda
H^{\psi}(Y^{f}(t))\right\\}=\mathsf{E}\exp\left\\{-\psi(\Lambda)Y^{f}(t)\right\\}=\tilde{\ell}_{f}(t,\psi(\Lambda)).$

To derive equation (36) we follow the same lines as is the proof of Theorem 2,
but we use now instead of the equation (22) (or (23)) for probabilites of the
process $N^{\psi}=N(H^{\psi}),$ the following equation for the probabilities
$\tilde{p}_{n}^{\psi}(t)$ of the process $M^{\psi}=M(H^{\psi})$(see [13]):

$\frac{d}{dt}\tilde{p}_{n}^{\psi}(t)=-\psi(\Lambda)\tilde{p}_{n}^{\psi}(t)-\sum_{m=1}^{n}\sum_{\Omega(k,m)}\psi^{z_{k}}(\Lambda)\prod_{j=1}^{k}\frac{(-\lambda_{j})^{x_{j}}}{x_{j}!}\tilde{p}_{n-m}^{\psi}(t).$

Thus, we obtain

${\mathbb{D}}_{t}^{f}\tilde{p}_{n}^{\psi,f}(t)=\psi(\Lambda)\tilde{p}_{n}^{\psi,f}(t)-\sum_{m=1}^{n}\sum_{\Omega(k,m)}\psi^{z_{k}}(\Lambda)\prod_{j=1}^{k}\frac{(-\lambda_{j})^{x_{j}}}{x_{j}!}\tilde{p}_{n-m}^{\psi,f}(t)+\tilde{p}_{n}^{\psi,f}(0)\nu(t),$

and then apply the same reasonings as those after formula (29) to come to
(36).

Using the expression for the probability generating function
$\tilde{G}^{\psi}(u,t)$ of the time-changed process
$M^{\psi}(t)=M(H^{\psi}(t))$ (see [11]):

$\tilde{G}^{\psi}(u,t)=\exp\Big{\\{}-t\psi\Big{(}\sum_{j=1}^{k}\lambda_{j}(1-u^{j})\Big{)}\Big{\\}},$

we calculate $\tilde{G}^{\psi,f}$ as follows:

$\displaystyle\tilde{G}^{\psi,f}(u,t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\tilde{G}^{\psi}(u,s)\ell_{f}(t,s)ds=$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\exp\Big{\\{}-t\psi\Big{(}\sum_{j=1}^{k}\lambda_{j}(1-u^{j})\Big{)}\Big{\\}}\ell_{f}(t,s)ds=$
$\displaystyle=$
$\displaystyle\mathsf{E}e^{-\psi\left(\sum_{j=1}^{k}\lambda_{j}(1-u^{j})\right)Y^{f}(t)}=\tilde{\ell}_{f}\Big{(}t,\psi\Big{(}\sum_{j=1}^{k}\lambda_{j}(1-u^{j})\Big{)}\Big{)},$

that is, (37) holds. Equation (38) follows in view of (37) and (9). ∎

######  Remark 3.

Following the same calculation as in [13, formula (5.3)], equation (36) can
also be written in the form

$\mathcal{D}^{f}_{t}\tilde{p}_{n}^{\psi,f}(t)=-\psi\Big{(}\Lambda\Big{(}I-\frac{1}{\Lambda}\sum_{j=1}^{n\wedge
k}\lambda_{j}B^{j}\Big{)}\Big{)}\tilde{p}_{n-m}^{\psi,f}(t).$

Example 2. As a continuation of Example 1, consider GCP with time-change given
by a compound Poisson-Gamma subordinator: $M^{GN}(t)=M(G_{N}(t))$.
Probabilities of this process can be calculated and expressed in terms of
special functions:

$\displaystyle\tilde{p}_{n}^{GN}(t)$ $\displaystyle=$
$\displaystyle\mathsf{P}\left\\{M(G_{N}(t))=n\right\\}=\int_{0}^{\infty}\mathsf{P}(M(s)=n)\mathsf{P}(G_{N}(t)\in
ds)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\left(\lambda_{j}s\right)^{x_{j}}}{x_{j}!}e^{-\Lambda
s}$ $\displaystyle\times\Big{\\{}e^{-\lambda
t}\delta_{\\{0\\}}(ds)+e^{-\lambda t-\beta
s}\frac{1}{s}\Phi\left(\alpha,0,\lambda t(\beta s)^{\alpha}\right)\Big{\\}}ds$
$\displaystyle=$ $\displaystyle e^{-\lambda
t}\int_{0}^{\infty}\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\left(\lambda_{j}s\right)^{x_{j}}}{x_{j}!}e^{-\Lambda
s}e^{-\beta s}\frac{1}{s}\Phi\left(\alpha,0,\lambda t(\beta
s)^{\alpha}\right)ds$ $\displaystyle=$ $\displaystyle e^{-\lambda
t}\sum_{\Omega(k,n)}\int_{0}^{\infty}\prod_{j=1}^{k}\frac{\left(\lambda_{j}s\right)^{x_{j}}}{x_{j}!}e^{-(\Lambda+\beta)s}\frac{1}{s}\sum_{l=1}^{\infty}\frac{(\lambda
t\beta^{\alpha})^{l}}{l!\Gamma(\alpha l)}s^{\alpha l-1}ds$ $\displaystyle=$
$\displaystyle e^{-\lambda
t}\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\sum_{l=1}^{\infty}\frac{(\lambda
t\beta^{\alpha})^{l}}{l!\Gamma(\alpha
l)}\int_{0}^{\infty}e^{-(\Lambda+\beta)s}s^{\sum_{j=1}^{k}x_{j}+\alpha l-1}ds$
$\displaystyle=$ $\displaystyle e^{-\lambda
t}\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\sum_{l=1}^{\infty}\frac{(\lambda
t\beta^{\alpha})^{l}}{l!\Gamma(\alpha l)}\frac{\Gamma\left(z_{k}+\alpha
l\right)}{\left(\Lambda+\beta\right)^{z_{k}+\alpha l}}$ $\displaystyle=$
$\displaystyle e^{-\lambda
t}\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\frac{1}{\left(\Lambda+\beta\right)^{z_{k}}}\sum_{l=1}^{\infty}\frac{(\lambda
t\beta^{\alpha})^{l}}{l!\left(\Lambda+\beta\right)^{\alpha
l}}\frac{\Gamma\left(z_{k}+\alpha l\right)}{\Gamma(\alpha l)},$

which can be written in the form

$\tilde{p}_{n}^{GN}(t)=e^{-\lambda
t}\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\frac{1}{\left(\Lambda+\beta\right)^{z_{k}}}\,{{}_{1}}\Psi_{1}\left((z_{k},\alpha),(0,\alpha),\frac{\lambda
t\beta^{\alpha}}{\left(\Lambda+\beta\right)^{\alpha}}\right),$

where ${{}_{p}}\Psi_{q}$ is the generalized Wright function

${}_{p}\Psi_{q}((a_{i},\alpha_{i}),(b_{j},\beta_{j}),z)=\sum_{k=0}^{\infty}\frac{\prod\limits_{i=1}^{p}\Gamma(a_{i}+\alpha_{i}k)}{\prod\limits_{j=1}^{q}\Gamma(b_{j}+\beta_{j}k)}\frac{z^{k}}{k!}$

defined for $z\in C$, $a_{i},b_{i}\in C$, $\alpha_{i},\beta_{i}\in R$,
$\alpha_{i},\beta_{i}\neq 0$ and $\sum\alpha_{i}-\sum\beta_{i}>-1$ (see, e.g.,
[9]).

If $\alpha=1$, that is, $G_{N}(t)=E_{N}(t)$, then probabilities become of the
form

$\tilde{p}_{n}^{E}(t)=e^{-\lambda
t}\sum_{\Omega(k,n)}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\frac{1}{\left(\Lambda+\beta\right)^{z_{k}}}\mathcal{E}_{1,2}^{z_{k}+1}\left(\frac{\lambda
t\beta}{\Lambda+\beta}\right),$

where $\mathcal{E}_{\rho,\delta}^{\gamma}$ is the three-parameter generalized
Mittag-Leffler function, which is defined as follows:

$\mathcal{E}_{\rho,\delta}^{\gamma}(z)=\sum_{k=0}^{\infty}\frac{\Gamma(\gamma+k)}{\Gamma(\gamma)}\frac{z^{k}}{k!\Gamma(\rho
k+\delta)},\text{ }z\in C,\text{ }\rho,\delta,\gamma\in C,$

with $\text{ Re}(\rho)>0,\text{Re}(\delta)>0,\text{Re}(\gamma)>0$ (see, e.g.,
[9]). For $n=0$ we have:

$\tilde{p}_{0}^{GN}(t)=e^{-\lambda t}+e^{-\lambda
t}\sum_{l=1}^{\infty}\frac{\left(\lambda
t\beta^{\alpha}\right)^{l}}{\left(\Lambda+\beta\right)^{\alpha
l}}\frac{\Gamma(\alpha l)}{l!\Gamma(\alpha l)}=\exp\Big{\\{}-\lambda
t\Big{(}1-\frac{\beta^{\alpha}}{\left(\Lambda+\beta\right)^{\alpha}}\Big{)}\Big{\\}},$

and if $\alpha=1$

$\tilde{p}_{0}^{E}(t)=\exp\Big{\\{}-\frac{\lambda\Lambda
t}{\Lambda+\beta}\Big{\\}}.$

Consider now the process $M^{GN,f}(t)=M(G_{N}(Y^{f}(t))$, $t\geq 0$. The next
theorem follows from Theorem 5.

######  Theorem 6.

The process $M^{GN,f}$ has probability distribution function given by formula
(35) with $\psi=\psi_{GN}$ and the probabilities $\tilde{p}_{n}^{GN,f}(t)$
satisfy the following equation:

$\displaystyle{\mathcal{D}}_{t}^{f}\tilde{p}_{n}^{GN,f}(t)$ $\displaystyle=$
$\displaystyle\left(\frac{\lambda\beta^{\alpha}}{\left(\Lambda+\beta\right)^{\alpha}}-\lambda\right)\tilde{p}_{n}^{GN,f}(t)+$
$\displaystyle+\frac{\left(\lambda\beta\right)^{\alpha}}{\left(\Lambda+\beta\right)^{\alpha}\Gamma(\alpha)}\sum_{m=1}^{n}\sum_{\Omega(k,m)}\frac{\Gamma\left(z_{k}+\alpha\right)}{\left(\Lambda+\beta\right)^{z_{k}}}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\tilde{p}_{n-m}^{GN,f}(t).$

The process $M^{E,f}$has probability distribution function given by (35) with
$\psi=\psi_{EN}$ and the probabilities $\tilde{p}_{n}^{E,f}(t)$ satisfy the
equation

$\displaystyle{\mathcal{D}}_{t}^{f}\tilde{p}_{n}^{E,f}(t)$ $\displaystyle=$
$\displaystyle-\lambda\frac{\Lambda}{\Lambda+\beta}\tilde{p}_{n}^{E,f}(t)+$
$\displaystyle+\frac{\lambda\beta}{\Lambda+\beta}\sum_{m=1}^{n}\sum_{\Omega(k,m)}\frac{\Gamma\left(z_{k}+1\right)}{\left(\Lambda+\beta\right)^{z_{k}}}\prod_{j=1}^{k}\frac{\lambda_{j}^{x_{j}}}{x_{j}!}\tilde{p}_{n-m}^{E,f}(t).$

## References

  * [1] Beghin, L., D’Ovidio, M. Fractional Poisson process with random drift. Electron. J. Probab. 19: 1–26 (2014)
  * [2] Beghin, L., Orsingher, E. Fractional Poisson processes and related planar random motions. Electron. J. Probab. 14, 1790–1826 (2009)
  * [3] Buchak, K., Sakhno, L. Compositions of Poisson and Gamma processes. Mod. Stoch.: Theory Appl. 4(2), 161–188 (2017)
  * [4] Buchak, K., Sakhno, L. Properties of Poisson processes directed by compound Poisson-Gamma subordinators. Mod. Stoch.: Theory Appl. 5(2), 167–189 (2018)
  * [5] Buchak, K., Sakhno, L. On the governing equations for Poisson and Skellam processes time-changed by inverse subordinators. Theor. Probab. Math. Stat. 98, 91–104 (2019)
  * [6] Chukova, S., Minkova, L. Pólya-Aeppli of order $k$ risk model. Comm. Statist. Simulation Comput. 44(3), 551–564 (2015)
  * [7] Di Crescenzo, A., Martinucci, B. and Meoli, A. A fractional counting process and its connection with the Poisson process. ALEA, Lat. Am. J. Probab. Math. Stat. 13, 291–307 (2016)
  * [8] Garra, R., Orsingher, E., Scavino, M. Some probabilistic properties of fractional point processes Stoch. Anal. Appl. 35(4), 701–718 (2017)
  * [9] Haubold, H.J., Mathai, A.M., Saxena, R.K. Mittag-Leffler functions and their applications. Journal of Applied Mathematics. 51 (2011)
  * [10] Kadankova, T., Leonenko, N. and Scalas, E. Fractional non-homogeneous Poisson and Pólya-Aeppli processes of order $k$ and beyond. Comm. Statist. Theory Methods 52(8), 2682–2701 (2023)
  * [11] Kataria, K.K., Khandakar M. Skellam and time-changed variants of the generalized fractional counting process. Fractional Calculus and Applied Analysis. 25, 1873–1907 (2022)
  * [12] Kataria, K.K., Khandakar M. Generalized Fractional Counting Process. Journal of Theoretical Probability 35(4), 2784–2805 (2022)
  * [13] Kataria, K.K., Khandakar M. Generalized Counting Process: its Non-Homogeneous and Time-Changed Versions. Preprint arXiv:2210.03981 (2022)
  * [14] Kochubei, A.N. General fractional calculus, evolution equations, and renewal processes. Integral Equ. Oper. Theory. 71, no. 4, 583–600 (2011)
  * [15] Kostadinova, K., Minkova, L. 2019. On the Poisson process of order $k$. Pliska Stud. Math. Bulgar. 22, 117–28 (2019)
  * [16] Leonenko, N., Scalas, E., Trinh, M. The fractional non-homogeneous Poisson process. Stat. Probab. Lett. 120, 147–156 (2017)
  * [17] Meerschaert, M.M., Toaldo, B. Relaxation patterns and semi-Markov dynamics. Stoch. Proc. Appl. 129, Issue 8, 2850–2879 (2019)
  * [18] Orsingher, E., Polito, F. The space-fractional Poisson process. Stat. Probab. Lett. 82, 852–858 (2012)
  * [19] Orsingher, E., Polito, F. On the integral of fractional Poisson processes. Stat. Probab. Lett. 83(4), 1006–1017 (2013)
  * [20] Orsingher, E., Toaldo, B. Counting processes with Bernštein intertimes and random jumps. J. Appl. Probab. 52, 1028–1044 (2015)
  * [21] Toaldo, B. Convolution-type derivatives, hitting-times of subordinators and time-changed $C_{0}$-semigroups. Potential Analysis. 42, 115–140 (2015)


11institutetext: College of Engineering, Mathematics & Physical Sciences,
University of Exeter, UK
22institutetext: Institute of Computer Science, Foundation for Research &
Technology - Hellas (FORTH), Greece

# Deep Semantic Manipulation of Facial Videos

Girish Kumar Solanki 11 0000-0002-4493-1726 Anastasios Roussos
2211 0000-0001-6015-3357<EMAIL_ADDRESS>
###### Abstract

Editing and manipulating facial features in videos is an interesting and
important field of research with a plethora of applications, ranging from
movie post-production and visual effects to realistic avatars for video games
and virtual assistants. Our method supports semantic video manipulation based
on neural rendering and 3D-based facial expression modelling. We focus on
interactive manipulation of the videos by altering and controlling the facial
expressions, achieving promising photorealistic results. The proposed method
is based on a disentangled representation and estimation of the 3D facial
shape and activity, providing the user with intuitive and easy-to-use control
of the facial expressions in the input video. We also introduce a user-
friendly, interactive AI tool that processes human-readable semantic labels
about the desired expression manipulations in specific parts of the input
video and synthesizes photorealistic manipulated videos. We achieve that by
mapping the emotion labels to points on the Valence-Arousal space (where
Valence quantifies how positive or negative is an emotion and Arousal
quantifies the power of the emotion activation), which in turn are mapped to
disentangled 3D facial expressions through an especially-designed and trained
expression decoder network. The paper presents detailed qualitative and
quantitative experiments, which demonstrate the effectiveness of our system
and the promising results it achieves.

##  1 Introduction

Manipulation and synthesis of photorealistic facial videos is a significant
challenge in computer vision and graphics. It has plenty of applications in
visual effects, movie post-production for the film industry, video games,
entertainment apps, visual dubbing, personalized 3D avatars, virtual reality,
telepresence and many more fields. Concerning photographs, commercial and
research software allows editing colors and tone of photographs [4, 18] and
even editing visual style [34]. Moreover, for manipulating facial videos, the
traditional methods involve the use of physical markers on the faces and
expensive setups of lights and cameras at different angles in controlled
conditions, along with the use of complex CGI and VFX tools [2]. More robust,
affordable, and data-driven approaches to learn from the facial features of
the subject would allow easy and fine-grained control over manipulations and
would be applicable even in videos captured “in-the-wild”.

Figure 1: Main steps of the proposed pipeline for semantic manipulation of
facial videos.

The vast majority of existing works perform manipulations of facial features
in images rather than videos and typically require manual input by an
experienced user. Some approaches are image-based and lack parametric control
[5, 19, 48], while other approaches like [42] perform a parametric control by
mapping parameters of a 3D Morphable Model (3DMM) to a generative model like
the so-called StyleGAN [23]. However, this disentangled control is on the
images generated by Generative Adversarial Networks (GANs) and not on the
original images. Tewari et al. [41] overcome this limitation by developing an
embedding for real portrait images in the latent space of the StyleGAN and
support photorealistic editing on real images. The so-called ICface [46] makes
use of Action Units (AUs) to manipulate expressions in images, producing good
results, which however are not as visually plausible as the results produced
by methods based on 3D models. Moreover, the manual control using AUs can be
cumbersome due to too many parameters and the requirement of expert knowledge.
Kollias et al. [27] synthesize visual affect from neutral facial images, based
on input in the form of either basic emotion labels or Valence-Arousal (VA)
pairs, supporting the generation of image sequences. Wang et al. [47]
synthesize talking-head videos from a single source image, based on a novel
keypoint representation and decomposition scheme. Even though this kind of
methods output photorealistic videos and produce promising results, they solve
a different problem than this paper, since they process input coming from just
a single image and do not attempt to perform multiple semantic manipulations
on a whole video of a subject. Gafni et al. [15] use neural radiance fields
for monocular 4D facial avatar reconstruction, which synthesize highly-
realistic images and videos, having the ability to control several scene
parameters. However, their method focuses on facial reenactment tasks and when
used to manipulate a subject’s facial video, the control over the facial
expressions is rather limited and relies again on altering the parameters of a
3DMM.

This paper proposes a user-driven and intuitive method that performs
photorealistic manipulation of facial expressions in videos. We introduce an
interactive AI tool for such manipulations that supports user-friendly
semantic control over expressions that is understandable to a wide variety of
non-experts. Our AI tool allows the user to specify the desired emotion label,
the intensity of expression, and duration of manipulation from the video
timeline and performs edits coherent to the entire face area. We propose a
novel robust mapping of basic emotions to expression parameters of an
expressive 3D face model, passing through an effective intermediate
representation in the Valence-Arousal (VA) space. In addition, we build upon
state-of-the-art neural rendering methods for face and head reenactment [30,
13], achieving photorealistic videos with manipulated expressions.

##  2 Related Work

####  2.0.1 Facial expression representation and estimation.

Thanks to the revolution that Deep Learning and Convolutional Neural Networks
(CNNs) have brought in the field, recent years have witnessed impressive
advancements in video-based facial expression analysis [35]. For example,
Barros et al. [6] use a deep learning model to learn the location of emotional
expressions in a cluttered scene. Otberdout et al. [36] use covariance
matrices to encode deep convolutional neural network (DCNN) features for
facial expression recognition and show that covariance descriptors computed on
DCNN features are more efficient than the standard classification with fully
connected layers and SoftMax. Koujan et al. [29] estimate 3D-based
representation of facial expressions invariant to other image parameters such
as shape and appearance variations due to identity, pose, occlusions,
illumination variations, etc. They utilize a network that learns to regress
expression parameters from 3DMMs and produces 28D expression parameters
exhibiting a wide range of invariance properties. Moreover, Toisoul et al.
[45] jointly estimate the expressions in continuous and categorical emotions.

3D Morphable Models (3DMMs) play a vital role in 3D facial modelling and
reconstruction. They are parametric models and capable of generating the 3D
representation of human faces. Following the seminal work of Blanz and Vetter
[7], the research in this field has been very active until today and during
the last decade 3DMMs have been successfully incorporated in deep learning
frameworks [14]. Some recent advancements in this field produced more powerful
and bigger statistical models of faces with rich demographics and variability
in facial deformations [22, 12, 31, 8].

Facial reenactment transfers facial expressions from a source to a target
subject to conditioning the generative process on the source’s underlying
video. Traditional methods of facial reenactment use either 2D wrapping
techniques [33, 16] or a 3D face model [44] manipulating only the face
interior. More recent works are typically based on GANs, which are deep
learning frameworks that have produced impressive results in image and video
synthesis [21]. Neural Textures [43] perform real-time animation of facial
expression of a target video by using an image interpolation technique to
modify deformations within the internal facial region of the target and
redirecting manipulated face back to the original target frame. There is no
control over the head pose and eye gaze in face reenactment systems, but they
are very useful in applications like video dubbing [17]. Kim et al. [24]
further improve the realism of video dubbing by preserving the style of a
target actor. They perform monocular reconstruction of source and target to
use expressions of source, preserving identity, pose, illumination, and eyes
of the target.

Many facial reenactment techniques exploit 3DMMs, since they offer a reliable
way to separate expressions and identity from each other. For example, Deep
Video Portraits [25] use a GAN-based framework to translate these 3D face
reconstructions to realistic frames, being the first to transfer the
expressions along with full 3D head position, head rotation, eye gaze and
blinking. However, they require long training times, and the results have an
unnatural look due to the inner mouth region and teeth. Head2Head [30]
overcomes these limitations by using a dedicated multiscale dynamics
discriminator to ensure temporal coherence, and a dedicated mouth
discriminator to improve the quality of mouth area and teeth. This work was
extended to Head2Head++[13] by improving upon the 3D reconstruction stage and
faster way to detect eye gaze based on 68 landmarks, achieving nearly real-
time performance. Papantoniou et al. [37] build upon Head2Head++ and introduce
a method for the photorealistic manipulation of the emotional state of actors
in videos. They introduce a deep domain translation framework to control the
manipulation through either basic emotion labels or expressive style coming
from a different person’s facial video. However, they do not exploit more
fine-grained and interpretable representations of facial expressions, such as
the VA representation adopted in our work.

Figure 2: GUI of our AI tool that processes semantic labels about the desired
expression manipulations in specific parts of the input video and synthesizes
photorealistic manipulated videos. The user can play the input and the
manipulated video side by side to observe the effect of the manipulations.
Please also refer to the supplementary video for a demonstration of the GUI
[1].

##  3 Proposed method

We introduce a system that estimates and manipulates disentangled components
of a 3D face from facial videos followed by self-reenactment, while
maintaining a high level of realism. As depicted in Fig. 1, for every specific
subject whose facial video we wish to manipulate, the main steps of our
pipeline are as follows (more details are given in the following subsections):

1) We collect several videos with a varied range of expressions for the
specific subject. 2) We process the collected videos and estimate for every
frame of every video a data pair of VA values and Expression coefficients of a
3D face model by applying EmoNet [45] and DenseFaceReg [13] respectively. 3)
We use the data pairs generated in step 2 to train a person-specific
Expression Decoder Network to learn a mapping from 2D VA values to 30D
expression coefficients. 4) For the video of the subject that we wish to
manipulate, we train a state-of-the-art head reenactment method [13] to
produce a neural rendering network that produces photorealistic videos of the
subject. 5) To support the determination of semantic labels by the user, we
assume input in the form of basic expression labels alongside their intensity,
which we map to the VA space and from there to 3D expression coefficients,
through the trained Expression Decoder Network (step 3). 6) We input the
manipulated 3D expression coefficients to the trained neural renderer (step
4), which outputs the manipulated video.

###  3.1 AI tool

As a Graphical User Interface (GUI) for the manipulations by the user, we
introduce an interactive AI tool that we call “Deep Facial Manipulator” (DFM),
Fig. 2. A user can open the video to be edited in the tool and can
interactively add the edits by selecting a semantic label for the emotion as
well as its intensity with start and end duration from the timeline. The
“Reenact Video” button from the tool runs the reenactment process and the
edited video is played in the tool once the process finishes. The user can
then add more edits or modify the existing ones and repeat the process. Once
the user is happy with the manipulations, they can keep the saved manipulated
video. Please refer to the supplementary video for a demonstration of the
usage of our AI tool [1]. The main steps that implement the AI tool’s
functionality are described in more detail in the following sections.

Figure 3: Construction of a Person-Specific Dataset of Facial Expressions. We
process a collection of videos with varied expressions and construct a large
set of data pairs, where each pair consists of a 2D VA vector and a 30D
expression coefficients vector.

###  3.2 Person-specific dataset of facial expressions

As already mentioned, for every subject whose video we want to manipulate, we
collect a set of additional facial videos of the same subject with a varied
range of expressions. This video collection helps our Expression Decoder
Network learn the way that the specific subject is expressing their emotions.
It is worth mentioning that the videos of this collection do not need to be
recorded under the same conditions and scene as the video we wish to
manipulate since we use and robustly estimate a disentangled representation of
3D facial expressions that is invariant to the scene conditions. For example,
in the case that the subject is a celebrity actor, we can collect short clips
from diverse interviews and different movie scenes with the specific actor. In
more detail, we process every frame of every video of the person-specific
collection to estimate the following (see Fig. 3):

Valence-arousal values: First, we perform face detection and facial
landmarking in every frame by applying the state-of-the-art method of [9].
Then, we resize the cropped face image to a resolution of 256x256 pixels and
feed it to the state-of-the-art EmoNet method [45] to estimate Valence-
Arousal values. We have chosen Valence-Arousal values as an intermediate
representation of facial expression because these provide a continuous
representation that matches more accurately the variability of real human
expressions, as compared to discrete emotion labels (like Happy, Sad, Angry,
etc.), while at the same time they can be easily associated to expression
labels specified by a nonexpert user of our AI tool. Fig. 4(a) shows an
example of the distribution of VA values for one of the videos of the
collection of a specific subject.

3D expression coefficients: For every frame, we perform 3D face reconstruction
by running the so-called DenseFaceReg network from [13]. This is a robust and
efficient CNN-based approach that estimates a dense 3D facial mesh for every
frame, consisting of about 5K vertices. Following [13], we fit to this mesh
the combined identity and expression 3DMM used in [30, 13] and keep the
30-dimensional expression coefficients vector that describes the 3D facial
deformation due to facial expressions in a disentangled representation. By
performing the aforementioned estimations in the person-specific video
collection, we construct a dataset of data pairs consisting of a 2D Valence-
Arousal vector and a 30D expression coefficients vector.

|
---|---
(a) | (b)
Figure 4: Adopted representation of facial expressions on the Valence-Arousal
space. (a) Exemplar scatter plot of the Valence-Arousal points from one of our
dataset videos. We observe that the points occupy a significant part of the VA
space. (b) Our pre-defined mapping from expression label and intensity to a
specific region of the VA space, which we use to convert the user input in our
AI tool. Please note that ”Neutral” is the only label that does not have
different strength levels.

###  3.3 Expression decoder network

Our Expression Decoder Network is a person-specific mapping from the Valence-
Arousal space to the expression coefficients of the 3D face model. For every
subject, we train it in a supervised manner using the person-specific dataset
generated in Section 3.2. In more detail, we create a multilayer perceptron
network that takes a 2D VA vector and regresses a 30D expression coefficient
vector at the output. The loss function minimized during training is the Root
Mean Square Error (RMSE) loss between the estimated and ground-truth
expression coefficients. The network consists of 6 fully-connected layers with
4096, 2048, 1024, 512, 128 and 64 units per layer respectively and with
Rectified Linear Units (RELU) to introduce non-linearities [20]. The number of
units and layers are decided based on the empirical testing on a smaller
sample of data with different configurations and was tuned by trying different
learning rates and batch sizes. A learning rate of $10^{-3}$ and batch size of
32 was selected, based on the trials. To avoid overfitting and generalize
better, “dropout” is used as a regularization technique where some network
units are randomly dropped while training, hence preventing units from co-
adapting too much [40]. Furthermore, we used Adam optimizer [26] for updating
gradients and trained the network for 1000 epochs. Figure 5 demonstrates an
example of results obtained by our Expression Decoder Network (please also
refer to the supplementary video [1]). The input VA values at different time
instances form a 2D orbit in the VA space. For visualization purposes, the
output expression coefficients provided by our Decoder Network are visualized
as 3D facial meshes by using them as parameters of the adopted 3DMM [30, 13]
in combination with a mean face identity. The visualized 3D meshes show that
our Decoder Network can generate a range of facial expressions in a plausible
manner.

Figure 5: Visualization of exemplar inference results of our Expression
Decoder Network. $1^{st}$ row: sequence of input VA values, that form a 2D
orbit, with the red dot signifying the current location. $2^{nd}$ row:
corresponding output expression coefficients, visualized as a 3D mean face
having the expression that is specified by these coefficients.

###  3.4 Synthesis of photorealistic manipulated videos

An overview of our system’s module for synthesizing photorealistic videos with
manipulated expressions is provided in Fig. 6. We are based on the framework
of GANs [21] and build upon the state-of-the-art method of Head2Head++ [13],
which produces high-quality results in head reenactment scenarios. We use the
video to be manipulated as training footage for Head2Head++ to train a neural
renderer that synthesizes controllable sequences of fake frames of the subject
in the video to be manipulated. The process involves the estimation of facial
landmarks and eye pupils, as well as disentangled components of the 3D face
(identity, expressions, pose), which helps us effectively modify the
expression component while keeping the other components unaltered. These
components are then combined to create a compact image-based representation
called Normalized Mean Face Coordinate (NMFC) image, which is further combined
with the eye gaze video and used as conditional input to the GAN-based video
renderer to synthesize the fake frames. Head2Head++ also takes care of
temporal coherence between frames and has a dedicated mouth discriminator for
a realistic reenactment of the mouth region.

Each facial expression manipulation is specified by the user in the user-
friendly format of a time interval of the edit, emotion label, and intensity
(low, medium, or high). Adopting one of the two sets supported by EmoNet [45],
we consider the following 5 basic emotion labels: neutral, happy, sad,
surprise, fear. Inspired by [45], we use a pre-defined mapping from emotion
label and intensity to a specific region of the VA space, see Fig. 4(b). To
further improve the realism of the results, we randomly sample VA values
within this specific region and associate them with equally spaced time
instances within the time interval. We use B-spline interpolation to connect
these VA values and create a sequence of VA values so that there is one VA
pair for every frame. These VA values are then fed to our Expression Decoder
Network to generate the manipulated expression parameters. To ensure a
perceptually plausible transition of expressions, in case that the desired
expression is different from neutral, we transit from neutral to the desired
expression and from the desired expression back to neutral using 20 frames at
the beginning and end of each edit. Furthermore, to also ensure smooth
dynamics of facial expressions and avoid any noise and jittering artifacts, we
use approximating cubic splines for smoothing [38].

Figure 6: Overview of our module for synthesis of photorealistic manipulated
videos.

Having computed a vector of expression coefficients for every frame to be
manipulated, we combine it with the identity and camera parameters from the
original frame to produce manipulated NMFC images, see Fig. 6. The sequence of
manipulated NMFC images is then utilized by the trained neural renderer to
synthesize photorealistic manipulated frames with the desired expressions.

Figure 7: Frames from exemplar manipulations using our deep semantic
manipulation system. The output manipulated videos correspond to several
different emotion labels.

##  4 Experimental results

We conduct qualitative as well as quantitative experiments to evaluate the
proposed system. The experiments presented in this section are conducted on
videos from 2 subjects: one male and one female (hereafter referred to as
subject A and subject B respectively - see Fig. 7). All used videos have a
frame rate of 30 fps. The video collection for the person-specific datasets of
facial expressions (Sec. 3.2) originated from videos shot with a mobile camera
for Subj. A (18 videos with an overall duration of 1 hour 40 mins) and from
publicly available YouTube videos for Subj. B (20 videos with an overall
duration of 1 hour 55 mins). These video collections included variability in
terms of facial expressions, to train our Expression Decoder Network.
Furthermore, for each subject, the adopted neural renderer was trained on a
single video (with 2 mins 17 sec of duration for subj. A and 1 min 23 sec of
duration for subj. B), after following the pre-processing pipeline of [13],
i.e. face detection, image cropping and resizing to a resolution of 256x256
pixels. Please refer to the supplementary video for additional results and
visualizations [1].

###  4.1 Qualitative evaluation

We apply our method to generate several types of expression manipulations in
videos, altering for example the emotion class from happy to sad, from neutral
to happy, from sad to fear, and so on. Fig. 7 illustrates examples of video
frames with such manipulations, using different types of manipulated
expressions. We observe that our method succeeds in drastically changing the
depicted facial expressions in a highly photorealistic manner. In several
cases, the altered facial expressions seem to correspond well to the target
emotion labels (e.g., in the case of the label “happy”), but this is not
always the case (e.g., in the label “sad”). This can be attributed to the
different levels of availability and diversity of training data for different
emotions.

Figure 8: Frames from exemplar manipulations using our system, where the
target expressions and intensities are controlled. Please also see the
supplementary video [1].

Fig. 8 presents examples of controlling the intensity of the target
expressions, which is also supported by the proposed method. It visualizes
specific frames of the manipulated videos, with different manipulated
expressions at different intensities (low, medium, and high). We observe that
for the label “happy”, the effect of varying the controlled intensities is
visually apparent, whereas for the label “sad” all three intensity levels
yield very similar results. Furthermore, in the visualized example, the target
label “happy” matches very well the perceived emotion by humans.

###  4.2 Quantitative evaluation

Following [30, 13, 25] we quantitatively evaluate the synthesized videos
through self-reenactment, since this is the only practical way to have ground
truth at the pixel level. In more detail, the videos of subjects A and B that
have been selected for training the person-specific video renderer undergo a
70%-30% train-test split (the first 70% of the video duration is used for
training and the rest for testing). The frames at the train set are used to
train the adopted neural renderer as normal, whereas the frames at the test
set are used to evaluate the self-reenactment and are considered as ground
truth. The following pre-processing steps are applied to the test set: i) 3D
face reconstruction to extract disentangled components (expression
coefficients, identity, pose) for every frame, and ii) extraction of values
for every frame using EmoNet [45]. Subsequently, we test 2 types of self-
reenactment: A) using “Ground Truth” expression coefficients, by which we mean
that we feed our method with the expression coefficients extracted through the
pre-processing step (i) described above. This type of self-reenactment
corresponds to the one used in previous works [30, 13, 25]. B) using
expression coefficients synthesized by our Expression Decoder Network, after
this is fed with valence-arousal values extracted through the pre-processing
step (ii) described above. This is a new type of self-reenactment which we
term as “emotion self-reenactment”. Following [13], we use the following
metrics: 1) Average Pixel Distance (APD), which is computed as the average
L2-distance of RGB values across all spatial locations, frames and videos,
between the ground truth and synthesized data. 2) Face-APD, which is like APD
with the only difference being that instead of all spatial locations, it
considers only the pixels within a facial mask, which is computed using the
NMFC image. 3) Mouth-APD, which is the same as Face-APD but the mouth
discriminator of Head2Head++ [13].

Fig. 9 visualizes the synthesized images and APD metrics in the form of
heatmaps for some example frames from the test set of subject A. We observe
that our method manages to achieve accurate emotion self-reenactment. This is
despite the fact that the synthesized expression of every frame comes from our
Expression Decoder Network, using as sole information the valence-arousal
values and without any knowledge about how these are instantiated in terms of
facial geometry and deformations at the specific ground truth frame. We also
see that the synthesized images as well as the values and spatial distribution
of the APD errors is very close to the case where GT expressions are being
used (“with GT expressions”), which solves the substantially less challenging
problem of simple self-reenactment. Table 1 presents the overall quantitative
evaluation in terms of APD metrics averaged over all pixels frames and both
subjects. We observe that, unsurprisingly, as we move from APD to Face-APD and
then to Mouth-APD, the error metrics increase, since we focus on more and more
challenging regions of the face. However, all APD metrics of our method are
relatively low and consistently very close to the case of simple self-
reenactment (“with GT expressions”).

Figure 9: Example frames from the self-reenactment evaluation. Our actual
model is compared with the simplified case where the GT expressions are used.
Despite solving a more difficult task, we achieve very similar quality, both
visually and in terms of APD measures. We also observe that the synthesized
frames are perceptually very similar to the source (ground truth). Please also
see the supplementary video [1]. Table 1: Quantitative evaluation of self-
reenactment with GT expressions and with expression decoder network (our
actual model). Average pixel distance (APD) values over all frames and over
both subjects tested are reported. The considered range of pixel values is
[0,255].

###  4.3 User studies

Realism. We designed a comprehensive web-based user study to evaluate the
realism of the manipulated videos, as perceived by human observers. Following
[30, 25], the questionnaire of the user study included both manipulated and
real videos, with a randomly shuffled order and without revealing to the
participants which videos were real. In addition, the participants were asked
to watch each video once and rate its realism as per their first perception on
a Likert-type scale of 1 to 5, where 1 means “absolutely fake” and 5 means
“absolutely real”. Scores with values 4 or 5 are considered as corresponding
to realistic videos. In total, 24short videos (with subjects A and B having an
equal share of 12 videos each) were included in the questionnaire and shown to
every participant of the user study. The manipulated videos were of two types:
A) videos with single manipulation (SM), which corresponds to a single
manipulated expression at the central part of the video (60% of video
duration) with a forward and backward transition from the real footage (real
to fake and then back to real footage), and B) videos with double manipulation
(DM), which included two manipulated expressions with similar forth and back
transitions from real footage. In total, 21 participants took part in this
user study.

Table 2 presents the results of the user study on realism. We observe that for
Subj. A, we achieve a realism score (percentage of “real”) that is almost the
same with the score for real videos (77% versus 79%). This reflects the high
realism of the results that our system is able to achieve. We also observe
that for Subj. B, the performance of our system is relatively lower, which
might be attributed to the fact that the training footage for this subject was
relatively shorter. However, even for this subject, our system succeeds in
synthesizing videos that are perceived as real more than 50% of the times,
which is still a promising result. In addition, we see that videos with single
manipulation (SM) achieved better realism scores as compared with videos with
double manipulation (DM). This is due to the fact that keeping a high level of
realism becomes more challenging when multiple transitions to fake expressions
are included in the manipulated video.

Table 2: Results of the user study about the realism of manipulated videos. 21
users participated in the study and the number of answers per score is
reported. ”Real” corresponds to the percentage of scores with values 4 or 5.

AI tool. We designed another web-based user study, where we asked participants
to try our AI tool “DFM”. The participants were provided with a small overview
of the tool and were asked to operate it via a screen share session. The
participants used the DFM tool to generate the manipulated videos and observe
the results. Based on their user experience, they answered some questions
regarding the tool and technology they used, based on a Likert-type scale from
1 to 5. The questions were related to the user experience in terms of usage,
functionality, recommendation to a friend, and how much they were impressed by
using such tool and technology. In total, 5 users took part in this study. The
results of the user study for our AI tool are presented in Table 3. We observe
that all users responded very positively (scores 4 or 5) to the questions
related to the ease of usage, the functionality, the chances of recommending
the tool to a friend, and the excitement of the technology used. These results
show the potential impact that our novel framework for photorealistic video
manipulation can have. On the other hand, three users rated the speed and
performance as 3 which is logical as users like faster results, and our
implementation has not yet reached real-time speeds. In addition, the answers
in the question regarding design and interaction show that some participants
believe that there is a scope of improvement in the GUI design for a better
user experience.

Table 3:  Results of the user study about our AI tool. The users rated
different aspects of the tool on a Likert-type scale from 1 to 5. 5 users
participated and the number of answers per question and score is reported.

##  5 Conclusion and future work

We proposed a novel approach to photorealistic manipulation of facial
expressions in videos. Our system gives the ability to the user to specify
semantically meaningful manipulations in certain parts of the input video in
terms of emotion labels and intensities. We achieve this by mapping the
emotion labels to valence-arousal (VA) values, which in turn are mapped to
disentangled 3D facial expressions, through our novel Expression Decoder
Network. In addition, we build upon the latest advances in neural rendering
for photorealistic head reenactment. We also developed an easy-to-use
interactive AI tool that integrates our system and gives the ability to a non-
technical common user to access this functionality. The extensive set of
experiments performed demonstrates various capabilities of our system and its
components. Moreover, the qualitative and quantitative evaluations along with
the results of the user studies provide evidence about the robustness and
effectiveness of our method. For future work, we aim to increase the speed of
the manipulations using the tool and eliminate the long training times of
videos to be edited on the neural renderer. This seems achievable with recent
advancements in works like Head2HeadFS [11] that use few-shot learning.
Furthermore, it might be interesting to develop a “universal” expression
decoder network that will be able to achieve realistic results without needing
to rely on such an extensive person-specific training.
Note on social impact. Apart from the various applications with positive
impact on society and our daily lives, this kind of methods raise concerns,
since they could be misused in a harmful manner, without the consent of the
depicted individuals [10]. We believe that scientists working in these fields
need to be aware of and seriously consider these risks and ethical issues.
Possible countermeasures include contributing in raising public awareness
about the capabilities of current technology and developing systems that
detect deepfakes [32, 3, 28, 39].
Acknowledgments. A.Roussos was supported by HFRI under the ‘$1^{st}$ Call for
HFRI Research Projects to support Faculty members and Researchers and the
procurement of high-cost research equipment’ Project I.C.Humans, Number 91.

## References

  * [1] https://youtu.be/VIGFHaa1aIA
  * [2] Abouaf, J.: Creating illusory realism through vfx. IEEE Computer Graphics and Applications 20(04), 4–5 (2000)
  * [3] Amerini, I., Caldelli, R.: Exploiting prediction error inconsistencies through lstm-based classifiers to detect deepfake videos. In: Proceedings of the 2020 ACM workshop on information hiding and multimedia security. pp. 97–102 (2020)
  * [4] Ashbrook, S.: Adobe [r] photoshop lightroom. PSA Journal 72(12), 12–13 (2006)
  * [5] Averbuch-Elor, H., Cohen-Or, D., Kopf, J., Cohen, M.F.: Bringing portraits to life. ACM Transactions on Graphics (ToG) 36(6), 1–13 (2017)
  * [6] Barros, P., Parisi, G.I., Weber, C., Wermter, S.: Emotion-modulated attention improves expression recognition: A deep learning model. Neurocomputing 253, 104–114 (2017)
  * [7] Blanz, V., Vetter, T.: A morphable model for the synthesis of 3d faces. In: Proceedings of the 26th annual conference on Computer graphics and interactive techniques. pp. 187–194. ACM Press/Addison-Wesley Publishing Co. (1999)
  * [8] Booth, J., Roussos, A., Ponniah, A., Dunaway, D., Zafeiriou, S.: Large scale 3d morphable models. International Journal of Computer Vision 126(2), 233–254 (2018)
  * [9] Bulat, A., Tzimiropoulos, G.: How far are we from solving the 2d & 3d face alignment problem?(and a dataset of 230,000 3d facial landmarks). In: Proceedings of the IEEE International Conference on Computer Vision. pp. 1021–1030 (2017)
  * [10] Chesney, B., Citron, D.: Deep fakes: A looming challenge for privacy, democracy, and national security. Calif. L. Rev. 107, 1753 (2019)
  * [11] Christos Doukas, M., Rami Koujan, M., Sharmanska, V., Zafeiriou, S.: Head2headfs: Video-based head reenactment with few-shot learning. arXiv e-prints pp. arXiv–2103 (2021)
  * [12] Dai, H., Pears, N., Smith, W.A., Duncan, C.: A 3d morphable model of craniofacial shape and texture variation. In: Proceedings of the IEEE international conference on computer vision. pp. 3085–3093 (2017)
  * [13] Doukas, M.C., Koujan, M.R., Sharmanska, V., Roussos, A., Zafeiriou, S.: Head2head++: Deep facial attributes re-targeting. IEEE Transactions on Biometrics, Behavior, and Identity Science 3(1), 31–43 (2021). https://doi.org/10.1109/TBIOM.2021.3049576
  * [14] Egger, B., Smith, W.A., Tewari, A., Wuhrer, S., Zollhoefer, M., Beeler, T., Bernard, F., Bolkart, T., Kortylewski, A., Romdhani, S., et al.: 3d morphable face models—past, present, and future. ACM Transactions on Graphics (TOG) 39(5), 1–38 (2020)
  * [15] Gafni, G., Thies, J., Zollöfer, M., Nießner, M.: Dynamic neural radiance fields for monocular 4d facial avatar reconstruction. In: IEEE/CVF Conf. on Computer Vision and Pattern Recognition (CVPR) (2021), https://justusthies.github.io/posts/nerface/
  * [16] Garrido, P., Valgaerts, L., Rehmsen, O., Thormahlen, T., Perez, P., Theobalt, C.: Automatic face reenactment. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 4217–4224 (2014)
  * [17] Garrido, P., Valgaerts, L., Sarmadi, H., Steiner, I., Varanasi, K., Perez, P., Theobalt, C.: Vdub: Modifying face video of actors for plausible visual alignment to a dubbed audio track. In: Computer graphics forum. vol. 34, pp. 193–204. Wiley Online Library (2015)
  * [18] Gatys, L.A., Ecker, A.S., Bethge, M.: Image style transfer using convolutional neural networks. In: 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). pp. 2414–2423 (2016). https://doi.org/10.1109/CVPR.2016.265
  * [19] Geng, J., Shao, T., Zheng, Y., Weng, Y., Zhou, K.: Warp-guided gans for single-photo facial animation. ACM Transactions on Graphics (ToG) 37(6), 1–12 (2018)
  * [20] Glorot, X., Bordes, A., Bengio, Y.: Deep sparse rectifier neural networks. In: Proceedings of the fourteenth international conference on artificial intelligence and statistics. pp. 315–323. JMLR Workshop and Conference Proceedings (2011)
  * [21] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., Bengio, Y.: Generative adversarial nets. Advances in neural information processing systems 27 (2014)
  * [22] Huber, P., Hu, G., Tena, R., Mortazavian, P., Koppen, P., Christmas, W.J., Ratsch, M., Kittler, J.: A multiresolution 3d morphable face model and fitting framework. In: Proceedings of the 11th international joint conference on computer vision, imaging and computer graphics theory and applications (2016)
  * [23] Karras, T., Laine, S., Aila, T.: A style-based generator architecture for generative adversarial networks. In: Proceedings of the IEEE/CVF conference on computer vision and pattern recognition. pp. 4401–4410 (2019)
  * [24] Kim, H., Elgharib, M., Zollhöfer, M., Seidel, H.P., Beeler, T., Richardt, C., Theobalt, C.: Neural style-preserving visual dubbing. ACM Trans. Graph. 38(6) (nov 2019). https://doi.org/10.1145/3355089.3356500, https://doi.org/10.1145/3355089.3356500
  * [25] Kim, H., Garrido, P., Tewari, A., Xu, W., Thies, J., Niessner, M., Pérez, P., Richardt, C., Zollhöfer, M., Theobalt, C.: Deep video portraits. ACM Trans. Graph. 37(4) (jul 2018). https://doi.org/10.1145/3197517.3201283, https://doi.org/10.1145/3197517.3201283
  * [26] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. iclr. 2015. arXiv preprint arXiv:1412.6980 9 (2015)
  * [27] Kollias, D., Cheng, S., Ververas, E., Kotsia, I., Zafeiriou, S.: Deep neural network augmentation: Generating faces for affect analysis. International Journal of Computer Vision 128(5), 1455–1484 (2020)
  * [28] Korshunov, P., Marcel, S.: Vulnerability assessment and detection of deepfake videos. In: 2019 International Conference on Biometrics (ICB). pp. 1–6. IEEE (2019)
  * [29] Koujan, M.R., Alharbawee, L., Giannakakis, G., Pugeault, N., Roussos, A.: Real-time facial expression recognition “in the wild” by disentangling 3d expression from identity. In: 2020 15th IEEE International Conference on Automatic Face and Gesture Recognition (FG 2020). pp. 24–31. IEEE (2020)
  * [30] Koujan, M.R., Doukas, M.C., Roussos, A., Zafeiriou, S.: Head2head: Video-based neural head synthesis. In: 2020 15th IEEE International Conference on Automatic Face and Gesture Recognition (FG 2020). pp. 16–23. IEEE (2020)
  * [31] Li, T., Bolkart, T., Black, M.J., Li, H., Romero, J.: Learning a model of facial shape and expression from 4D scans. ACM Transactions on Graphics 36(6), 194:1–194:17 (Nov 2017), two first authors contributed equally
  * [32] de Lima, O., Franklin, S., Basu, S., Karwoski, B., George, A.: Deepfake detection using spatiotemporal convolutional networks. arXiv preprint arXiv:2006.14749 (2020)
  * [33] Liu, Z., Shan, Y., Zhang, Z.: Expressive expression mapping with ratio images. In: Proceedings of the 28th annual conference on Computer graphics and interactive techniques. pp. 271–276 (2001)
  * [34] Luan, F., Paris, S., Shechtman, E., Bala, K.: Deep photo style transfer. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 4990–4998 (2017)
  * [35] Martínez, B., Valstar, M.F., Jiang, B., Pantic, M.: Automatic analysis of facial actions: A survey. IEEE Transactions on Affective Computing 10, 325–347 (2019)
  * [36] Otberdout, N., Kacem, A., Daoudi, M., Ballihi, L., Berretti, S.: Deep covariance descriptors for facial expression recognition. In: BMVC (2018)
  * [37] Papantoniou, F.P., Filntisis, P.P., Maragos, P., Roussos, A.: Neural emotion director: Speech-preserving semantic control of facial expressions in” in-the-wild” videos. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 18781–18790 (2022)
  * [38] Pollock, D.: Smoothing with cubic splines. Handbook of Time Series Analysis, Signal Processing, and Dynamics pp. 293–332 (12 1999). https://doi.org/10.1016/B978-012560990-6/50013-0
  * [39] Rössler, A., Cozzolino, D., Verdoliva, L., Riess, C., Thies, J., Nießner, M.: Faceforensics++: Learning to detect manipulated facial images. 2019 IEEE/CVF International Conference on Computer Vision (ICCV) pp. 1–11 (2019)
  * [40] Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., Salakhutdinov, R.: Dropout: a simple way to prevent neural networks from overfitting. The journal of machine learning research 15(1), 1929–1958 (2014)
  * [41] Tewari, A., Elgharib, M., Bernard, F., Seidel, H.P., Pérez, P., Zollhöfer, M., Theobalt, C.: Pie: Portrait image embedding for semantic control. ACM Transactions on Graphics (TOG) 39(6), 1–14 (2020)
  * [42] Tewari, A., Elgharib, M., Bharaj, G., Bernard, F., Seidel, H.P., Pérez, P., Zollhofer, M., Theobalt, C.: Stylerig: Rigging stylegan for 3d control over portrait images. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 6142–6151 (2020)
  * [43] Thies, J., Zollhöfer, M., Nießner, M.: Deferred neural rendering: Image synthesis using neural textures. ACM Transactions on Graphics (TOG) 38(4), 1–12 (2019)
  * [44] Thies, J., Zollhofer, M., Stamminger, M., Theobalt, C., Nießner, M.: Face2face: Real-time face capture and reenactment of rgb videos. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 2387–2395 (2016)
  * [45] Toisoul, A., Kossaifi, J., Bulat, A., Tzimiropoulos, G., Pantic, M.: Estimation of continuous valence and arousal levels from faces in naturalistic conditions. Nature Machine Intelligence 3(1), 42–50 (2021)
  * [46] Tripathy, S., Kannala, J., Rahtu, E.: Icface: Interpretable and controllable face reenactment using gans. In: Proceedings of the IEEE/CVF winter conference on applications of computer vision. pp. 3385–3394 (2020)
  * [47] Wang, T.C., Mallya, A., Liu, M.Y.: One-shot free-view neural talking-head synthesis for video conferencing. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (2021)
  * [48] Zakharov, E., Shysheya, A., Burkov, E., Lempitsky, V.: Few-shot adversarial learning of realistic neural talking head models. In: Proceedings of the IEEE/CVF international conference on computer vision. pp. 9459–9468 (2019)